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Chapter 1 

Introduction 



An outline for this course. 

• We will observe that many phenomena in ecology, biology and biochemistry can be 
modelled mathematically. 

• We will initially focus on systems where the spatial variation is not present or, at 
least, not important. Therefore only the temporal evolution needs to be captured 
in equations and this typically (but not exclusively) leads to difference equations 
and/or ordinary differential equations. 

• We are inevitably confronted with systems of non-linear difference or ordinary dif- 
ferential equations, and thus we will study analytical techniques for extracting in- 
formation from such equations. 

• We will proceed to consider systems where there is explicit spatial variation. Then 
models of the system must additionally incorporate spatial effects. 

• In ecological and biological contexts the main physical phenomenon governing the 
spatial variation is typically, but again not exclusively, diffusion. Thus we are in- 
variably required to consider parabolic partial differential equations. Mathematical 
techniques will be developed to study such systems. 

• These studies will be in the context of ecological, biological and biochemical appli- 
cations. In particular we will draw examples from: 

— enzyme dynamics and other biochemical reactions; 

— epidemics; 

— interaction ecological populations, such as predator-prey models; 

— biological pattern formation mechanisms; 

— chemotaxis; 

— the propagation of an advantageous gene through a population; 

— nerve pulses and their propagation. 
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1.1 References 

The main references for this lecture course will be: 

• J. D. Murray, Mathematical Biology, 3rd edition, Volume I [8]. 

• J. D. Murray, Mathematical Biology, 3rd edition, Volume II [9]. 
Other useful references include (but are no means compulsory): 

• J. P. Keener and J. Sneyd, Mathematical Physiology [7]. 

• L. Edelstein-Keshet, Mathematical Models in Biology [2]. 

• N. F. Britton, Essential Mathematical Biology [1]. 



Chapter 2 

Spatially independent models for a 
single species 

In this chapter we consider modelhng a single species in cases where spatial variation is not 
present or is not important. In this case we can simply examine the temporal evolution 
of the system. 

References. 

• J. D. Murray, Mathematical Biology, 3rd edition. Volume I, Chapter 1 and Chapter 
2 [8]. 

• L. Edelstein-Keshet, Mathematical Models in Biology, Chapter 1, Chapter 2 and 
Chapter 6 [2]. 

• N. F. Britton, Essential Mathematical Biology, Chapter 1 [1]. 

2.1 Continuous population models for single species 

A core feature of population dynamics models is the conservation of population number, 
i.e. 

rate of increase of population = birth rate — death rate (2-1) 

+ rate of immigration — rate of emigration. 

We will make the assumption the system is closed and thus there is no immigration or 
emigration. 

Let N{t) denote the population at time t. Equation (2.1) becomes 

^ = fiN) = Ng{N), (2.2) 
where g{N) is defined to be the intrinsic growth rate. Examples include: 
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The Malthus model. This model can be written as: 



g{N) = b-d''^^r, (2.3) 
where b and d are constant birth and death rates. Thus 

^ = rN, (2.4) 

and hence 

N{t) = Noe^'. (2.5) 
The Verhulst model. This model is also known as the logistic growth model: 

giN) =r(l-^Y (2.6) 



Definition. In the logistic growth equation, r is defined to be the linear birth rate 
and K is defined to be the carrying capacity. 



For N <^ K, we have 

~ riV ^ AT ~ ATqc^*. (2.7) 

However, as N tends towards K, 

dN , , 

-^0, (2.8) 

the growth rate tends to zero. 
We have 

dN ( N\ , , 

(2.9) 

and hence 

NnKe''^ 

N{t) = — -f— -^K as t^oo. (2.10) 

K + Noie^-^-l) ^ ' 

Sketching N{t) against time yields solution as plotted in Figure 2.1: we see that solutions 
always monotonically relaxes to as t — t- oo. 

Aside. The logistic growth model has been observed to give very good fits to popula- 
tion data in numerous, disparate, scenarios ranging from bacteria and yeast to rats and 
sheep [8]. 

2.1.1 Investigating the dynamics 

There are two techniques we can use to investigate the model 

^ = f{N) = Ng{N). (2.11) 
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Figure 2.1: Logistic growtli for Nq < K (left-liand) and Nq > K (riglit-liand). Parameters 
are as follows: r = 0.015 and K ~ 100. 



Method (i): analytical solution 

For the initial conditions N{t = 0) = A^O) with A'o fixed, we can we formally integrate 
equation (2.2) to give N(t) = N*{t), where N*{-) is the inverse of the function F{-) defined 
by 



F{x) 



1 



■ds. 



(2.12) 



However, unless integrating and finding the inverse function is straightforward, there is an 
easier way to determine the dynamics of the system. 

Method (ii): plot the graph 

Plot dN/dt = f{N) = Ng{N) as a function of N. For example, with 

f{N) = Ng{N) = N{6N'^ - - UN + 6) = N{N - 1){N - 2)(3 - N), (2.13) 
we have the plot shown in Figure 2.2. 



t 




1 2 
population (N) 



Figure 2.2: Growth according to the dynamics f{N) = N{N - 1){N - 2) (3 - N). 
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Note 1. For a given initial condition, Nq, the system will tend to the nearest root of 
f{N) = Ng{N) in the direction of /{Nq). The value of |iV(t)| will tend to infinity with 
large time if no such root exists. 

For f{N) = Ng{N) = N{N - 1){N - 2)(3 - N), we have: 

• when A'^o ^ (0,2] the large time asymptote is N{oo) = 1; 

• for A'^o > 2 the large time asymptote is N(oo) = 3; 

• N{t) = Vt if N{0) = 0. 



Note 2. On more than one occasion we will have a choice between using a graphical 
method and an analytical method, as seen above. The most appropriate method to use 
is highly dependent on context. The graphical method, Method (ii), quickly and simply 
gives the large time behaviour of the system and stability information (see below). The 
analytical method, Method (i), is often significantly more cumbersome, but yields all 
information, at a detailed quantitative level, about the system. 

Definition. A stationary point, also known as an equilibrium point, is a point where 
the dynamics does not change in time. Thus in our specific context of dN/dt = f{N) = 
Ng{N), the stationary points are the roots of f{N) = 0. 



Example. For dN/dt = f{N) = Ng{N) = N{N - 1){N - 2)(3 - A^), the stationary 
points are 

A^ = 0,1,2,3. (2.14) 



Definition. A stationary point is stable if a solution starting sufficiently close to the 
stationary point remains close to the stationary point. 



Non-examinable. A rigorous definition is as follows. Let A^Ary(t) denote the solution 
to dN/dt = f{N) = Ng{N) with initial condition N{t = 0) = A^q. A stationary point, 
A'^s, is stable if, and only if, for all e > there exists a 6 such that if [A'^s — A'ol < 5 then 
\NN,{t)-Ns\ <e. 



Exercise. Use Figure 2.2 to deduce which of stationary points of the system 
dN 



f{N) = Ng{N) = N{N - 1){N - 2)(3 - A^), (2.15) 



are stable. 



Solution. Figure 2.2 shows that both A^^ = 1 and A^^ = 3 are stable. 



Chapter 2. Spatially independent models for a single species 



11 



2.1.2 Linearising about a stationary point 

Suppose Ns is a stationary point of dN/dt = f{N) and make a small perturbation about 

N{t) = Ns + 7i{t), n{t)<^Ns. (2.16) 
We have, by using a Taylor expansion of f{N) and denoting ' = d/dN, that 

f{N{t)) = f{N, + n{t)) = f{N,) + n{t)f'{N,) + ^n{tff"{Ns) + ..., (2.17) 
and hence 

^ = ^ = f{N{t)) = f{Ns) + n{t)f'{N,) + ln{t)'f"{Ns) + ... (2.18) 

The linearisation of dN/dt = f{N) about the stationary point Ns is given by neglecting 
higher order (and thus smaller) terms to give 

^=m)n(,). 



The solution to this linear system is simply 



n{t) = n(t = 0) exp 



t^iNs) 
dN^ 



(2.19) 



Definition. Let Ns denote a stationary point of dN/dt = f{N), and let 



n{t) = n{t = 0) exp 



(2.20) 



be the solution of the linearisation about Ns- Then Ns is linearly stable if n{t) — ?• as 
t — 7- oo. In other words, Ns is linearly stable if 

^(A^.)<0. (2.21) 



Exercise. By algebraic means, deduce which stationary points of the system 
dN 

— = f{N) = Ng{N) = N{N - l)(iV - 2)(3 - N), (2.22) 
are linearly stable. Can your answer be deduced graphically? 
Solution. Differentiating f{N) with respect to gives 

f{N) = 2- 22N + ISiV^ - 4iV^ (2.23) 
and hence /'(O) = 6 (unstable), /'(I) = —8 (stable) etc. 

Consider the graph of f{N) to deduce stability graphically — steady states with negative 
gradient are linearly stable c.f. Figure 2.2. 
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population (N) 

Figure 2.3: Growth according to the dynamics f{N) = (1 — 7V)^. 



Exercise. Find a function f{N) such that dN/dt = f{N) has a stationary point which 
is stable and not linearly stable. 



Solution. The function 

f{N) = {l-Nf, 
gives /'(I) = and is therefore not linearly stable (see Figure 2.3). 



(2.24) 



2.1.3 Insect outbreak model 

First introduced by Ludwig in 1978, the model supposes budworm population dynamics 
to be modelled by the following equation: 



(2.25) 



The function p{N) is taken to represent the effect upon the population of predation by 
birds. Plotting p{N) as a function of N gives the graph shown in Figure 2.4. 




250 
population (N) 



500 



Figure 2.4: Predation, p{N), in the insect outbreak model. Parameters are as follows: 
A = 150, B = 0.5. 
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Non-dimensionsionalisation 

Let 



N = N*u, t = Tt, (2.26) 



where A^*, N have units of biomass, and t, T have units of time, with N* , T constant. 
Then 

N*du ( N*u\ B{N*fu^ 

-Yd^ = "^^H ^^J"^4^T(ivW ^ ^ ^ 

du ( N*u\ BTN*u'^ 

^ -r = rsTu 1 - -— - TTT-^. 2.28 

dr V Kb J + {N*Yv? ^ ' 

Hence with 

AT* A rj. ^ rp rsA Kb Kb , . 

N =A T=-, r = rBT = ^, ^ = _ = (2.29) 

we have 

du f 1 ^e/ 



Thus we have reduced the number of parameters in our model from four to two, which 
substantially simplifies our subsequent study. 

Steady states 

The steady states are given by the solutions of 

\ 2 
U\ U 



ru\\--\- -— ^ = 0. (2.31) 
V q) 1 + 

Clearly w = is a steady state. We proceed graphically to consider the other steady states 
which are given by the intersection of the graphs 

/i(^z)=r(l--') and f^{u) = (2.32) 

The top left plot of Figure 2.5 shows plots of fi{u) and /2('u) for different values of r and 
q. We see that, depending on the values of r and q, we have either one or three non-zero 
steady states. Noting that 

d/(u;r, q) 



= r>0, (2.33) 

u=0 



du 

typical plots of du/dr vs. u are shown in Figure 2.5 for a range of values of r and q 



Definition. A system displaying hysteresis exhibits a response to the increase of a 
driving variable which is not precisely reversed as the driving variable is decreased. 



Remark. Hysteresis is remarkably common. Examples include ferromagnetism and 
elasticity, amongst others. See http://en.wikipedia.org/wiki/Hysteresis for more 
details. 
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Figure 2.5: Dynamics of the non-dimensional insect outbreak model. Top left: plots of 
the functions fi{u) (dashed line) and f2{u) (solid line) with parameters r = 0.2,0.4,0.6, 
q = 10, 15, 20, respectively. Top right: plot of f{u; r, q) with parameters r = 0.6, q = 0.6. 
Bottom left: plot of /(w; r, q) with parameters r = 0.6, q = Q. Bottom right: plot of /(u; r, q) 
with parameters r = 0.6, q = 10. 



Extended Exercise 

• Fix r = 0.6. Explain how the large time asymptote of u, and hence N, changes as 
one slowly increases q from q <^ 1 to q ^ 1 and then one decreases q from g ^ 1 to 
g <C 1. In particular, show that hysteresis is present. Note for this value of r, there 
are three non-zero stationary points for q £ {qi, q2) with 1 < qi < q2 < 10. 

Solution. For small values of q there is only one non-zero steady state, S\. As q is 
increased past qi, three non-zero steady states exist. Si, S2, S3, but the system stays 
at Si. As q is increased further, past q2, the upper steady state ^3 is all that remains 
and hence the system moves to 53. If g is now decreased past 52, three non-zero 
steady states (5*1, ^2, ^3) exist but the system remains at ^3 until q is decreased 
past qi. 

Figure 2.6 shows f{u; r, q) for different values of q. The dashed line shows a plot for 
q = qi whilst the dash-dotted line shows a plot for q = q2- 

• What is the biological interpretation of the presence of hysteresis in this model? 
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Figure 2.6: Left-hand plot: dw/dr = /(u; r, q) in the non-dimensional insect outbreak model 
as q is varied. For small q there is one, small, steady state, for q G (91,92) there are three 
non-zero steady states and for large q there is one, large, steady state. Right-hand plot: the 
steady states plotted as a function of the parameter q reveals the hysteresis loop. 



Solution. If the carrying capacity, q, is accidentally manipulated such that an out- 
break occurs (Si —7- 5*3) then reversing this change is not sufficient to reverse the 
outbreak. 

2.1.4 Harvesting a single natural population 

We wish to consider a simple model for the maximum sustainable yield. Suppose, in the 
absence of harvesting, we have 

We consider a perturbation from the non-zero steady state, N = K. Thus we write 
= K + and find, on linearising, 

— - = -rn =^ n = noe"""*. (2.35) 

Hence the system returns to equilibrium on a timescale of 7^(0) = 0{l/r). 

We consider two cases for harvesting: 

• constant yield, Y; 

• constant effort, E. 

Constant yield 

For a constant yield, Y = Yq, our equations are 

!!^ = riv(l-^)-r„^l//(A';r„). (2,36) 
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25 50 75 100 

population (N) 

Figure 2.7 : Dynamics of the constant yield model for Yq = 0.00, 0.15, 0.30. As Yq is increased 
beyond a critical value the steady states disappear and — > in finite time. Parameters are 
as follows: K = 100 and r = 0.01. 



Plotting dN/dt as a function of reveals (see Figure 2.7) that the steady states disappear 
as Iq is increased beyond a critical value, and then A^ — )• in finite time. 
The steady states are given by the solutions of 

rN*-^-Yo = N* = ^ ^^^^ °^ . (2.37) 

Therefore extinction will occur once 

rK 

Yo > — . (2.38) 

Constant effort 

For harvesting at constant effort our equations are 

dN ( N\ def rN"^ 

^ = rN[l--yEN f{N; E) = N{r - E) - (2.39) 

where the yield is Y{E) = EN. The question is: how do we maximise Y{E) such that the 
stationary state still recovers? 

The steady states, N* , are such that f{N*;E) = (see Figure 2.8). Thus 

N*{E)=^^^:^^=(i-^)k, (2.40) 

and hence 

Y*{E) = EN*{E) = (^1- KE. (2.41) 

Thus the maximum yield, and corresponding value of A^*, are given by the value of E such 
that 

dY* r rK K 

-^ = ^ E = -, Y:,,, = —, N*^,, = -. (2.42) 
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Linearising about the stationary state N*{E) we have = N*{E) + n with 

n + . . . = — (r — E)n + . . . , 



dn , df{N;E) 

— ~ fp N ) H — — ■ — - 

dt ~ ' ^ dN 

and hence the recovery time is given by 



N=N* 



Tr{E) ~ o 



r-E 



(2.43) 



(2.44) 



Defining the recovery time to be the time for a perturbation to decrease by a factor of e 
according to the hnearised equations about the non-zero steady state, then 



rfl(o) = 

Hence, at the maximum yield state, 



1 



Tr{E) 



1 



r-E 



2 r 
Tr(E) = — since E = — at maximum yield. 
r 2 



(2.45) 



(2.46) 



As we measure Y it is useful to rewrite E in terms of Y to give the ratio of recovery times 
in terms of the yield Y{E) and the maximum yield Ym- 



Tb{Y) 
Tr{0) 

Derivation. At steady state, we have 
K 



(2.47) 



E^ - KE + Y* = as Y* = EN* =KE{1 
r \ r 



E 



(2.48) 



This gives 



E 



±r^l-AY*/K7 



r-E 



2 2 
Substituting into equation (2.45) gives the required result 



itWi 



(2.49) 



Plotting Tji(Y)/Tr{0) as a function of Y/Ym yields some interesting observations, as 
shown in Figure 2.8. 

Note. As Tr increases the population recovers less quickly, and therefore spends more 
time away from the steady state, N* . The biological implication is that, in order to 
maintain a constant yield, E must be increased. This, in turn, implies Tr increases, 
resulting in a positive feedback loop that can have disastrous consequences upon the 
population. 
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Figure 2.8: Dynamics of the constant effort model. The Icft-liand plot shows the logistic 
growth curve (solid line) and the yield, Y = EN (dashed lines), for two values of E. The 
right-hand plot shows the ratio of recovery times, Tii{Y)/Tii{0), with the negative root plotted 
as a dashed line and the positive root as a solid line. Parameters are as follows: K = 100 and 
r = 0.01. 



2.2 Discrete population models for a single species 

When there is no overlap in population numbers between each generation, we have a 
discrete model: 

Nt+i = Ntf{Nt) = H{Nt). (2.50) 

A simple example is 

Nt+i = rNt, (2.51) 

which implies 

oo r > 1 

Nt = r*iVo ^ "! No r = l . (2.52) 
r < 1 



Definition. An equilibrium point, N*, for a discrete population model satisfies 

N* = N*f{N*) = H{N*). (2.53) 
Such a point is often known fixed point. 



An extension of the simple model, equation (2.51), called the Ricker model includes a 
reduction of the growth rate for large Nt: 

Nt' ' 



Nt+i = Nt exp 
or, in non-dimensionalised form, 



r 1 



K 



, r > K>0, 



def 



ut+i=utex-p[r{l-ut)] = H{ut). 



(2.54) 
(2.55) 



We can start developing an idea of how this system evolves in time via cobwebbing, a 
graphical technique, as shown in Figure 2.9. 
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Figure 2.9: Dynamics of tlie Ricker model. The left-hand plot shows a plot of Nt+i = 
Nt exp [r (1 — Nt/K)] alongside iVt+i = Nt with the cobwebbing technique shown. The right- 
hand plot shows Nt for successive generation times t = 1, 2, . . . , 10. Parameters are as follows: 
iVo = 5, r = 1.5 and K = 100. 

In particular, it is clear that the behaviour sufficiently close to a fixed point, u* , depends 
on the value of H'{u*). For example: 

• -1< H'{u*) < 




t 



• H'{u*) = -1 




t 
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• H'{u*) < -1 




t 



2.2.1 Linear stability 

More generally, to consider the stability of an equilibrium point algebraically, rather than 
graphically, we write 

ut = u* + vt, (2.56) 

where u* is an equilibrium value. Note that u* is time-independent and satisfies u* = 
H{u*). Hence 

ut+i = u* + vt+i = H{u* + Vt) = H{u*) + vtH\u*) + o{v^). (2.57) 



Consequently, we have 

vt+i = H'{u*)vt where H'{u*) is a constant, independent of t, (2.58) 

and thus 

vt=[H'{u*)\\Q. (2.59) 
This in turn enforces stability if \H'{u*)\ < 1 and instability if \H'{u*)\ > 1. 

Definition. A discrete population model is linearly stable if \H'{u*)\ < 1. 
2.2.2 Further investigation 

The equations are not as simple as they seem. For example, from what we have seen thus 
far, the discrete time logistic model seems innocuous enough. 

Nt+i =rNt(^l-^^, r > K > 0. (2.60) 

If we put in enough effort, one could be forgiven for thinking that the use of cobwebbing 
will give a simple representation of solutions of this equation. However, the effects of 
increasing r are stunning. Figure 2.10 shows examples of cobwebbing when r = 1.5 and 
r = 4.0. 

It should now be clear that even this simple equation does not always yield a simple 
solution! How do we investigate such a complicated system in more detail? 
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Definition. A bifurcation point is, in the current context, a point in parameter space 
where the number of equilibrium points, or their stability properties, or both, change. 

We proceed to take a closer look at the non-dimensional discrete logistic growth model: 

ut+i = rut{l - ut) = H{ut), (2.61) 

for different values of the parameter r, and, in particular, we seek the values where the 
number or stability nature of the equilibrium points change. Note that we have equilibrium 
points at M* = and u* = (r — l)/r, and that H'{u) = r — 2ru. 

For < r < 1, we have: 

• = is a stable steady state since |-ff'(0)| = \r\ < 1; 

• the equilibrium point at u* = (r — l)/r is unstable. It is also unreachable, and thus 
irrelevant, for physical initial conditions with no > 0. 

For 1 < r < 3 we have: 

• u* = is an unstable steady state since |ff'(0)| = |r| > 1; 

• u* = {r — l)/r is an stable steady state since |if'((r — l)/r)| = |2 — r| < 1. 

In Figure 2.11 we plot this on a diagram of steady states, as a function of r, with stable 
steady states indicated by solid lines and unstable steady states by dashed lines. 

When r = 1 we have (r — l)/r = 0, so both equilibrium points are at u* = 0, with 
H'[u* = 0) = 1. Clearly we have a switch in the stability properties of the equilibrium 
points, and thus r = 1 is a bifurcation point. 
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Figure 2.11: Bifurcation diagram for the non-dimensional discrete logistic model. The non- 
zero steady state is given, for r > 1, by N* = {r — l)/r. 



What happens for r > 3? We have equihbrium points at n* = 0, u = (r — l)/r and 
H'{u* = (r — l)/r) < —1; both equilibrium points are unstable. Hence if the system 
approaches one of these equilibrium points the approach is only transient; it quickly moves 
away. We have a switch in the stability properties of the equilibrium points, and thus r = 3 
is a bifurcation point. 

To consider the dynamics of this system once r > 3 we consider 

Ut+2 = H{ut+i) = H [H{ut)] =^ H2{ut) = r [rut{l - Uf)] [1 - rut{l - ut)] . (2.62) 

Figure 2.12 shows H2{ut) for r = 2.5 and r = 3.5 and demonstrates the additional steady 
states that arise as r is increased past r = 3. 




Figure 2.12: Dynamics of the non-dimensional discrete logistic model in terms of every 
second iteration. The left-hand plot shows results for r = 2.5 whilst the right-hand plot shows 
results for r = 3.5. 



Note. The fixed points of H2 satisfy U2 = H2{u2), which is a quartic equation in 
However, we know two solutions, the fixed points H{-), i.e. and (r — l)/r. Using standard 
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techniques we can reduce the quartic to a quadratic, which can be solved to reveal the 
further fixed points of H2, namely 

ul = '^±^[{r-lf-A]"\ (2.63) 
Zr Zr ^ J 

These roots exist if (r — 1)^ > 4, i.e. r > 3. 



Definition. The m}^ composition of the function H is given by 



Hm{u)'^= [H ■ H ...H ■ H]{u). (2.64) 

^ V ' 

m times 





Definition. A point u is periodic of period 
Hm{u) = u, Hi{u) / tt, 


m for the function H if 
i € {1,2, ... m - 1}. 


(2.65) 



Thus the points 

u^ = !:±l±l[(r-l)^-4]^/^ (2.66) 
are points of period 2 for the function H. 

Problem. Show that the U2 are stable with respect to the function H2 for r > 3, 
(r-3) < 1. 

Let 

^^dgr + l^l 2_ 1,2^ ni = F(uo), U2 = H2{u,), (2.67) 

Zr Zr ^ J 

and let 



d_ 

du 

Then 



\ = — [H2{u)]\u=u,. (2.68) 



\ = ^[H- H{u)] \u=uo = H'{uo)H'{ui). (2.69) 

Thus for stability we require \H'{uq)H'{ui)\ < 1. 
Exercise. Finish the problem: show that the steady states 



are stable for the dynamical system ut+i = H2{ut), with r > 3, (r — 3) <C 1. 
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Exercise. Suppose uq is an equilibrium point of period m for the function H. Show 
that uq is stable if 

ni^o^ [H'{u,)] < 1, (2.71) 
where Ui = Hi{uo) for i G {1, 2, . . . , m — 1}. 



Solution. Defining A in a similar manner as before, we have 

d 



A 



du 
d_ 
du 



u=uo 

[H{Q{u)] 



where Q{u) = Hm-i{u), 



U=Uq 



d 

= H (itm-l)^-f^m-l 

Hence, by iteration, we have the result. 



(2.72) 
(2.73) 
(2.74) 
(2.75) 



u=uo 



We plot the fixed points of H2, which we now know to be stable, in addition to the fixed 
points of Hi, in Figure 2.13. The upper branch, U2u, is given by the positive root of 
equation (2.70) whilst the lower branch, U21, is given by the negative root. We have 
= H{u2u), u^u = H{u2i). Thus a stable, period 2, oscillation is present, at least for 
(r — 3) ^ 1. Any solution which gets sufficiently close to either U2u or ut^^^ stays close. 




Figure 2.13: Bifurcation diagram for the non-dimensional discrete logistic model with inclu- 
sion of the period 2 solutions. 



For higher values of r, there is a bifurcation point for H2; we can then find a stable 
fixed point for H/i[u) : : H2[H2{u)] in a similar manner. Increasing r further there is a 
bifurcation point for H4^(u). Again, we are encountering a level of complexity which is too 
much to deal with our current method. 

To bring further understanding to this complex system, we note the following definition. 
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Definition. An orbit generated by the point uq are the points {uo,ui,U2, ,U3, . . .} 
where Ui = Hi{u) = H{ui-i). 

We are primarily interested in the large time behaviour of these systems in the context of 
biological applications. Thus, for a fixed value of r, we start with a reasonable initial seed, 
say u* = 0.5, and plot the large time asymptote of the orbit of m*, ie. the points Hi{u*) 
once i is sufficiently large for there to be no transients. This gives an intriguing plot; see 
Figure 2.14. 

1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

3.0 3.2 3.4 3.6 3.8 4.0 
r 

Figure 2.14: The orbit diagram of the logistic map. For each value of r € [3,4] along 
horizontal axis, points on the large time orbits of the logistic map arc plotted. 




In particular, we have regions where, for r fixed, there are three points along the ver- 
tical corresponding to period 3 oscillations. This means any period of oscillation exists 
and we have a chaotic system. This can be proved using Sharkovskii's theorem. See P. 
Glendinning, Stability, Instability and Chaos [4] for more details on chaos and chaotic 
systems. 

Note. A common discrete population model in mathematical biology is 

Nt+i = * , . (2.76) 

Models of this form for the Colorado beetle are within the periodic regimes, while Nichol- 
son's blowfly model is in the chaotic regime [8]. 

2.2.3 The wider context 

In investigating the system 

ut+i = rut (1 - ut) H{ut), (2.77) 

we have explored a very simple equation which, in general, exhibits greatly different be- 
haviours with only a small change in initial conditions or parameters {i.e. linear growth 
rate, r). Such sensitivity is a hallmark feature of chaotic dynamics. In particular, it makes 
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prediction very difficult. There will always be errors in a model's formalism, initial condi- 
tions and parameters and, in general, there is no readily discernible pattern in the way the 
model behaves. Thus, assuming the real system behaviour is also chaotic, using statistical 
techniques to extract a pattern of behaviour to thus enable an extrapolation to predict 
future behaviour is also fraught with difficulty. Attempting to make accurate predictions 
with models containing chaos is an active area of research, as is developing techniques to 
analyse seemingly random data to see if such data can be explained by a simple chaotic 
dynamical system. 



Chapter 3 



Continuous population models: 
interacting species 



In this chapter we consider interacting populations, but again in the case where spatial 
variation is not important. Appendix A contains relevant information for phase plane 
analysis that may be useful. 

References. 

• J. D. Murray, Mathematical Biology, 3rd edition, Volume I, Chapter 3 [8]. 

• L. Edelstein-Keshet, Mathematical Models in Biology, Chapter 6 [2]. 

• N. F. Britton, Essential Mathematical Biology, Chapter 2 [1]. 

There are three main forms of interaction: 

Predator-prey An upsurge in population I (prey) induces a growth in population II 
(predator). An upsurge in population II (predator) induces a decline in population 
I (prey). 

Competition An upsurge in either population induces a decline in the other. 

Symbiosis An upsurge in either population induces an increase in the other. 

Of course, there are other possible interactions, such as cannibalism, especially with the 
"adult" of a species preying on the young, and parasitism. 

3.1 Predator-prey models 

The most common predator-prey model is the Lotka-Volterra model. With the number 
of prey and P the number of predators, this model can be written 



dA^ 

dP 
dF 



cNP - dP, 



aN - bNP, 



(3.1) 



(3.2) 
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with a, b, c, d positive parameters and c < b. Non-dimensionalising with u = {c/d)N, 
V = {b/a)P, T = at and a = d/a, we have 

1 ddu ad bda du ^ rr x /o on 

— ; — = U —UV =^ —— = U — UV = U(l — V) = t(U,V), (o.o) 

l/acdT c c b dr ^ j j \ ^ \ j 

1 adv da a dv 

— ■j- — — = c——uv — d—v, =^ — = a\uv — V) = aviu — 1) = q{u,v), 3.4 
l/abdT cb b dr ^ ' ^ i ^ \ j 

There are stationary points at {u,v) = (0,0) and {u,v) = (1, 1). 

Exercise. Find the stabihty of the stationary points {u,v) = (0,0) and {u,v) = (1,1). 
The Jacobian, J, is given by 



fu fv \ ^ I 'i^-v -u 
Qu Qv I \ av a{u - 1) 



(3.5) 



At (0, 0) we have 



-a 

with eigenvalues 1, —a. Therefore the steady state (0,0) is an unstable saddle. 
At (1, 1) we have 



-1 
a 



(3.7) 

with eigenvalues ziz^^/a. Therefore the steady state (1, 1) is a centre (not linearly stable). 



These equations are special; we can integrate them once, as follows, to find a conserved 
constant: 

du u(l — v) /"u— 1, f 1 — V , . 

A- = -r — rf ^ / / — • 

dv a{u — l)v J u J av 

Hence 

H = const = au + V — alnu — Inv. (3.9) 

This can be rewritten as 

/ pV\ / pU\a^ 

' ^ ^ e^, (3.10) 



V / \ u 

from which we can rapidly deduce that the trajectories in the (n, v) plane take the form 
shown in Figure 3.1. Thus u and v oscillate in time, though not in phase, and hence we 
have a prediction; predators and prey population numbers oscillate out of phase. There 
are often observations of this e.g. hare-lynx interactions. 
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0.0 1.0 2.0 3.0 4.0 5.0 4 8 12 16 

u Time 

Figure 3.1: Dynamics of tlie non-dimensional Lotka-Volterra system for a = 1.095 and 
H — 2.1,2.4,3.0,4.0. The left-hand plot shows the dynamics in the {u,v) phase plane whilst 
the right-hand plot shows the temporal evolution of u and v. 



3.1.1 Finite predation 

The common predator-prey model assumes that as — ?• oo the rate of predation per 
predator becomes unbounded, as does the rate of increase of the predator's population. 
However, with an abundance of food, these quantities will saturate rather than become 
unbounded. Thus, a more realistic incorporation of an abundance of prey requires a refine- 
ment of the Lotka-Volterra model. A suitable, simple, model for predator-prey interactions 
under such circumstances would be (after a non-dimensionalisation) 

du „, s. , . auv 

— = f{u,v) = u{l-u)--—, 3.11 
dr a + u 

g{u,v)=bv(l--), (3.12) 



dr \ uy 

where a, b, d are positive constants. Note the effect of predation per predator saturates 
at high levels of u whereas the predator levels are finite at large levels of prey and drop 
exceedingly rapidly in the absence of prey. 

There is one non-trivial equilibrium point, (u*,?;*), satisfying 

OjU* 

V* = V* where il — u*) = , (3.13) 

and hence 

* 

u = - 
2 

The Jacobian at [u* ,v* ) is 



-(a + d - 1) + V(a + d - 1)2 + 4d . (3.14) 



( fu 


fu \ 


\ gu 


gv J 



(3.15) 



{u* ,v*) 



where 



fu{u ,v ) = l-2u - — — - + =-u + — — -2 . (3.16 

d + u* [d + u*)'^ [d + u*)'' 
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b(v*)^ 

9u{u*y) = ^^ = b, (3.18) 



g,{u\v*) = b[l-2^—j =-h. (3.19) 

The eigenvalues satisfy 

(A-/n)(A-5.)- A'5« = X^-{fu + 9v)X + {fu9v-fvgu) = 0, (3.20) 
and hence 



A^-aA + /3 = =, A= "^^f (3.21) 



where 



Note that 



0=1- «(^*)' = 1 _ = (n*+d)-n* + (u*)^ ^ d+(n*)^ > .3 

(n* + d)2 {u*+d) u*+d d + u* ■ > 

Thus, if a < the eigenvalues A are such that we have either: 

• a stable node (a^ - 4/3 > 0); 

• stable focus (q^ - 4/3 < 0); 
at the equilibrium point {u*,v*). 

If a > we have an unstable equilibrium point at (n* , f * ) . 

3.2 A look at global behaviour 

This previous section illustrated local dynamics: we have conditions for when the dynamics 
will stably remain close to the non-trivial equilibrium point. One is also often interested 
in the global dynamics. However, determining the global dynamics of a system, away 
from its equilibrium points, is a much more difficult problem compared to ascertaining the 
local dynamics, sufficiently close to the equilibrium points. For specific parameter values, 
one can readily solve the ordinary differential equations to consider the behaviour of the 
system. One is also interested in the general properties of the global behaviour. This is 
more difficult, and we will consider one possible approach below. 

There are many potential tools available: nullcline analysis, the Poincare-Bendixson The- 
orem, the Poincare Index and the Bendixson-Dulac Criterion. The Poincare-Bendixson 
Theorem is a useful tool for proving that limit cycles must exist, while Poincare indices 
and the Bendixson-Dulac Criterion are useful tools for proving a limit cycle cannot exist. 
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We will briefly consider nullclines and the Poincare-Bendixson Theorem in detail. Please 
refer to P. Glendinning, Stability, Instability and Chaos: An Introduction to the Theory 
of Nonlinear Differential Equations [4], or D. W. Jordan and P. Smith, Mathematical 
Techniques: An Introduction for Engineering, Mathematical and Physical Sciences [6], for 
further details than considered here. 

3.2.1 Nullclines 



Definition. Consider the equations 

^ = f{u,v), f^=9{n,v). (3.24) 
The nullclines are the curves in the phase plane where f{u, v) = and g{u, v) = 0. 



Reconsider 



du X / N auu 

— = f{u,u) = u{l-u)- -—, 3.25 
dr a + u 

dv 
d7 



g{u,v) =bv(l-^y (3.26) 



The u nullclines are given by 



f(u,v) = u = and v = -(1 - u)(u + d). (3.27) 

a 

The V nullclines are given by 

g{u,v) = v = and v = u. (3.28) 

A sketch of the nullclines and the behaviour of the phase plane trajectories is shown in 
Figure 3.2. 

3.2.2 The Poincare-Bendixson Theorem 

For a system of two first order ordinary differential equations, consider a closed bounded 
region D. Suppose a positive half path, H, lies entirely within D. Then one of the 
following is true: 

1. H is a closed trajectory, e.g. a limit cycle; 

2. H asymptotically tends to a closed trajectory, e.g. a limit cycle; 

3. H terminates on a stationary point. 

Therefore, if D does not have a stationary point then there must be a limit cycle. 

For a proof see P. Glendinning, Stability, Instability and Chaos: An Introduction to the 
Theory of Nonlinear Differential Equations [4] . 
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Figure 3.2: The (u, v) phase-plane for the finite predation model when the steady state is 
stable. The u nullclines are plotted in red and the v nuUclines in green. Trajectories for a 
number of different initial conditions are shown as dashed lines. Parameters are as follows: 
a = 2.0, b = 0.1, d 2.0. 

Exercise. Explain why a > in the previous example (see equation (3.22)) implies 
we have limit cycle dynamics. What does this mean in terms of the population levels of 
predator and prey? 

Solution. For a > the steady state is an unstable node or spiral. Further, we can find a 
simple, closed boundary curve, C, in the positive quadrant of the {u,v) plane, such that 
on C phase trajectories always point into the domain, T>, enclosed by C. Applying the 
Poincare-Benedixon Theorem to the domain gives the existence of a limit cycle. See J. D. 
Murray, Mathematical Biology Volume I [8] (Chapter 3.4) for more details. 

3.3 Competitive exclusion 

We consider an ordinary differential equation model of two competitors. An example might 
be populations of red squirrels and grey squirrels [8]. Here, both populations compete for 
the same resources and a typical model for their dynamics is 



and N2 with grey squirrels in our example. 

In particular, given a range of parameter values and some initial values for A'^i and at 
the time t = 0, we would typically like to know if the final outcome is one of the following 
possibilities: 




(3.29) 



(3.30) 



where Ki, K2, ri, r2, 612, ^21 are positive constants. Let us associate A''i with red squirrels 



• the reds become extinct, leaving the greys; 
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• the greys become extinct, leaving the reds; 

• both reds and greys become extinct; 

• the reds and greys co-exist. If this system is perturbed in any way will the reds and 
greys continue to coexist? 

After a non-dimensionalisation (exercise) we have 

u'l = ui{l-ui-ai2U2)'^= fi{ui,U2), (3.31) 
U2 = PU2{1-U2- a2lUi)'^= f2{ui,U2), (3.32) 

where p = r2/ri. 

The stationary states are 

(nt,n^) = (0,0), «,^) = (1,0), (nt,^) = (0,1), (3.33) 

and ^ 

{ul,u*2) = - (1 - ai2, 1 - 021), (3.34) 

1 — ai2Q;2i 

if ai2 < 1 and 021 < 1 or Q12 > 1 and a2i > 1- 
The Jacobian is 



1 - 2ui — Qi2ti2 — Qi2ttl 

-pa2iU2 p{l - 2u2 - a2iui) 



(3.35) 



It is a straightforward application of phase plane techniques to investigate the nature of 
these equilibrium points: 



Steady state (nj,n2) = (0,0). 

(' 

p- X 



J - AI = ( ^ ^ ^ \ X = l,p. (3.36) 



Therefore (0, 0) is an unstable node. 
Steady state (nj,n2) = (1,0). 

J-AI=("\"^ , ,) A = -l,Kl-a2i). (3.37) 

y p{i- a2i) -Ay 

Therefore (1,0) is a stable node if a2i > 1 and a saddle point if a2i < 1. 
Steady state (n^,n2) = (0, 1). 

J - AI = I ^ " " ^ ° I A = -p, 1 - ai2. (3.38) 

\^ -yoa2i -p - X J 

Therefore (0, 1) is a stable node if ai2 > 1 and a saddle point if ai2 < 1. 
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Steady state K,n^) = Tr^7^(l - ai2,l - 021] 



J-AI 



1 / "21 — 1 — A ai2(ai2 — 1) 
1 - 012^21 y pa2i(a2i - 1) p{a2i - 1) - A 

Stability depends on Q12 and Q21. 



(3.39) 



There are several different possible behaviours. The totality of all behaviours of the above 
model is reflected in how one can arrange the nullclines within the positive quadrant. 
However, for competing populations these straight line nullclines have negative gradients. 
Figure 3.3 shows the model behaviour for different sets of parameter values. 




0.0 0.4 



1.2 



0.8 



0.4 



0.0 



/ 


/ / 




/ 


/ / 




^ / 


/ ' 

/ y 






















/ 






/ 












^ " " \Nv 


— 



0.0 



0.4 



1.2 




Figure 3.3: Dynamics of the non-dimensional competitive exclusion system. Top left: Q12 = 
0.8 < 1, a2\ — 1.2 > 1 and is excluded. Top right: a\2 — 1.2 > 1, 0,21 — 0.8 < 1 and u\ is 
excluded. Bottom left: ayi = 1.2 > 1, ol2\ = 1-2 > 1 and exclusion is dependent on the initial 
conditions. Bottom right: oi\2 = 0.8 < 1, ol2\ = 0.8 < 1 and we have coexistence. The stable 
steady states are marked with *'s and p = 1.0 in all cases. The red lines indicate /i = whilst 
the green lines indicate /2 = 0. 



Note. In ecology the concept of competitive exclusion is that two species competing for 
exactly the same resources cannot stably coexist. One of the two competitors will always 
have an ever so slight advantage over the other that leads to extinction of the second 
competitor in the long run (or evolution into distinct ecological niches). 
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3.4 Mutualism (symbiosis) 



We consider the same ordinary differential equation model for two competitors, i.e. 

dNi 



diV2 / N2 ^, iVi' 



(3.40) 
(3.41) 



where Ki, K2, ri, r2, 612, &21 are positive constants or, after non-dimensionalisation, 

def 



Ui{l - Ui + ai2U2) = fliui,U2), 
def 



PU2{1 - U2 + a2lUi] 



f2{ui,U2). 



(3.42) 
(3.43) 



In symbiosis, the straight line nullclines will have positive gradients leading to the following 
two possible behaviours shown in Figure 3.4. 





Figure 3.4: Dynamics of the non-dimensional symbiotic system. The left-hand figure shows 
population explosion (ai2 = 0.6 = 0:21) whilst the right-hand figure shows population coex- 
istence (ai2 = 0.1 — a2i). The stable steady states are marked with *'s and p = 1.0 in all 
cases. The red lines indicate /i = whilst the green lines indicate /2 = 0. 



3.5 Interacting discrete models 

It is also possible, and sometimes useful, to consider interacting discrete models which 
take the form 

ut+i = f{ut,vt), (3.44) 
vt+i = g{ut,vt), (3.45) 

and possess steady states at the solutions of 

u* = f{u*,v*), v* = g{u*,v*). (3.46) 

It is interesting and relevant to study the linear stability of these equilibrium points, and 
the global dynamics, but we do not have time to pursue this here. 



Chapter 4 

Enzyme kinetics 



In this chapter we consider enzyme kinetics, which can be thought of as a particular case 
of an interacting species model. In all cases here we will neglect spatial variation. 

Throughout, we will consider the m chemical species Ci, . . . , Cm- 

• The concentration of Cj, denoted q, is defined to be the number of molecules of Ci 
per unit volume. 

• A standard unit of concentration is moles m~^, often abbreviated to mol m~^. Recall 
that 1 mole = 6.023 x 10^'^ molecules. 

References. 

• J. D. Murray, Mathematical Biology, 3rd edition, Volume I, Chapter 6 [8]. 

• J. P. Keener and J. Sneyd, Mathematical Physiology, Chapter 1 [7]. 

4.1 The Law of Mass Action 

Suppose Ci, . . . , Cm undergo the reaction 

XlCi + X2C2 + . . . + XmCm V UiCi + U2C2 + . . . + l^mCm- (4.1) 

h 

The Law of Mass Action states that the forward reaction proceeds at rate 

kj4^4\..ct, (4.2) 

while the back reaction proceeds at the rate 

kf,C^^ ■ ■ ■ Cm ^ (4-3) 

where kj and kh are dimensional constants that must be determined empirically. 
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Note 1. Strictly, to treat kj, kb above as constant, we have to assume that the tem- 
perature is constant. This is a very good approximation for most biochemical reactions 
occurring in, for example, physiological systems. However, if one wanted to model re- 
actions that produce extensive heat for example, burning petrol, one must include the 
temperature dependence in /cj and kf, and subsequently keep track of how hot the sys- 
tem gets as the reaction proceeds. This generally makes the modelling significantly more 
difficult. Below we assume that we are dealing with systems where the temperature is 
approximately constant as the reaction proceeds. 

Note 2. The Law of Mass Action for chemical reactions can be derived from statistical 
mechanics under quite general conditions (see for example L. E. Riechl, A Modern Course 
in Statistical Physics [11]). 

Note 3. As we will see later, the Law of Mass Action is also used in biological scenarios to 
write down equations describing, for example, the interactions of people infected with, and 
people susceptible to, a pathogen during an epidemic. However, in such circumstances its 
validity must be taken as an assumption of the modelling; in such scenarios one cannot rely 
on thermodynamic/statistical mechanical arguments to justify the Law of Mass Action. 



4.2 Michaelis-Menten kinetics 

Michaelis-Menten kinetics approximately describe the dynamics of a number of enzyme 
systems. The reactions are 

ki 

S + E . SE, (4.4) 

k-i 

k2 

SE ^ P + E. (4.5) 

Letting c denoting the concentration of the complex SE, and s, e, p denoting the con- 
centrations of S, E, P, respectively, we have, from the Law of Mass Action, the following 
ordinary differential equations: 

ds 

— = —kise + k^ic, (4.6) 
dc 

— = kise — k^ic — k2C, (4.7) 
de 

— = —kise + k-ic + k2C, (4-8) 
I = k,c. (4.9) 

Note that the equation for p decouples and hence we can neglect it initially. 
The initial conditions are: 



s{0) = so, e(0) = eo < sq, c(0) = 0, p(0) = 0. (4.10) 
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Key Point. In systems described by the Law of Mass Action, linear combinations of 
the variables are often conserved. In this example we have 

A(e + c) = ^ e = eo-c, (4.11) 
at 

and hence the equations simplify to: 

ds 

— = -ki{eo - c)s + k_ic, (4.12) 
dc 

— = ki{eo - c)s - {k_i + k2)c, (4.13) 
with the determination of p readily achievable once we have the dynamics of s and c. 

4.2.1 Non-dimensionalisation 

We non-dimensionalise as follows: 

, , S C k2 def eo def k^^rk2 , . 

T = kieot, u = — , V = — , A = , e = — < 1, K = — ; , (4.14) 

So eo fciSQ So kiso 



which yields 



u' = -u + {u + K - X)v, (4.15) 
ev' = u-{u + K)v, (4.16) 



where u(0) = 1, v{0) = and e ^ 1. Normally e ~ 10 ^. Setting e = yields 

u 



(4.17) 



u + K' 

which is inconsistent with the initial conditions. Thus we have a singular perturbation 
problem; there must be a (boundary) region with respect to the time variable around t = 
where v' oo 0{\). Indeed for the initial conditions given we find v'(0) ~ 0(l/e), with m(0), 
f (0) < 0(1). This gives us the scaling we need for a singular perturbation investigation. 

4.2.2 Singular perturbation investigation 

We consider 

a = ^, (4.18) 

with 

u(t, e) = u{(T,e) = uo{(7) + eui{(T) + . . . , (4.19) 

v{t,€) = v{a,e)=VQ{a) + evi{a) + .... (4.20) 

Proceeding in the usual way, we find that uq, vq satisfy 

— = =^ uq = constant = 1, (4.21) 
dcr 

and 

^=uo-{l + K)vo = l-{l + K)vo vo = -— , (4.22) 

dcr 1 + K 
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which gives us the inner solution. 

To find the outer solution we expand 

u(r,e) = 
v{T,e) = 



within the equations 



to find that 



and 



This gives 



no(T) + eni(T) + . . . , 
Voir) + evi{T) + . . . , 



ev 

dug 

dr 



= -u+{u + K-X)v, 
= u — {u + K)v, 

= -no + {uq + K - X)vo, 
= no - (no + i^)no. 



no dno Ano 

— and — — = 

uq + K dr 



no + i^ 

In order to match the solutions as cr — )• oo and r — )• we require 



lim no = hm no = 1 and lim no = 1™ vq = — . 

(T— ^OO T— >0 IT— >0O T— >-0 1 -\- K 



(4.23) 
(4.24) 



(4.25) 
(4.26) 

(4.27) 
(4.28) 

(4.29) 



(4.30) 



Thus the solution looks like that shown in Figure 4.1. 




0.1 0.2 
time (t) 
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\- \ 
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0.5 1.0 
time (t) 



1.5 



Figure 4.1: Numerical solution of the non-dimensional Michaelis-Menten equations clearly 
illustrating the two different time scales. The u dynamics arc indicated by the solid line and 
the V dynamics by the dashed line. Parameters are e = 0.01, K = 0.1 and A = 1.0. 
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Often the initial, fast, transient is not seen or modelled and one considers just the outer 
equations with a suitably adjusted initial condition (ultimately determined from consis- 
tency/matching with the inner solution). Thus one often uses Michaelis-Menten kinetics 
where the equations are simply: 

du Xu , . N , u , . 

— = with u(0) = 1 and v = -. (4.31) 

dt u + K ^ ^ u + K ^ ^ 



Definition. We have, approximately, that dv/dr ~ using Michaelis-Menten kinet- 
ics. Taking the temporal dynamics to be trivial, 

^ 0, (4.32) 

when the time derivative is fast, i.e. of the form 

e^=giu,v), (4.33) 

where e ^ 1, g{u,v) ~ 0{1), is known as the pseudo-steady state hypothesis and is a 
common assumption in the literature. We have seen its validity in the case of enzyme 
kinetics about at least away from the inner region. 



Note. One must remember that the Michaelis-Menten kinetics derived above are a very 
useful approximation, but that they hinge on the validity of the Law of Mass Action. 
Even in simple biological systems the Law of Mass Action may breakdown. One (of 
many) reasons, and one that is potentially relevant at the sub-cellular level, is that the 
system in question has too few reactant molecules to justify the statistical mechanical 
assumptions underlying the Lass of Mass Action. Another reason is that the reactants are 
not well-mixed, but vary spatially as well as temporally. We will see what happens in this 
case later in the course. 



4.3 More complex systems 

Here we consider a number of other simple systems involving enzymatic reactions. In 
each case the Law of Mass Action is used to write down a system of ordinary differential 
equations describing the dynamics of the various reactants. See J. Keener and J. Sneyd, 
Mathematical Physiology [7], for more details. 

4.3.1 Several enzyme reactions and the pseudo-steady state hypothesis 

We can have multiple enzymes. In general the system of equations reduces to 



= f{u,Vi,...,Vn), 
= gi{u,Vi,...,Vn), 



(4.34) 
(4.35) 
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for i G {1, . . . , n}, while the pseudo-steady state hypothesis gives a single ordinary differ- 
ential equation 

u' = f{u,vi{u), . . . ,Vn{u)), (4.36) 
where vi{u), . . . , Vn{u) are the appropriate roots of the equations 

gi{u,vi, . . . ,Vn) = 0, i£{l,...,n}. (4.37) 

4.3.2 Allosteric enzymes 

Here the binding of one substrate molecule at one site affects the binding of another 
substrate molecules at other sites. A typical reaction scheme is: 

S + E . Ci > P + E (4.38) 

S + Ci ^— ^ C2 > Ci + E. (4.39) 

k-3 

Further details on the investigation of such systems can be found in J. D. Murray, Mathe- 
matical Biology Volume I [8], and J. P. Keener and J. Sneyd, Mathematical Physiology [7]. 

4.3.3 Autocatalysis and activator-inhibitor systems 

Here a molecule catalyses its own production. The simplest example is the reaction scheme 

A + b\2B, (4.40) 

though of course the positive feedback in autocatalysis is usually ameliorated by inhibition 
from another molecule. This leads to an example of an activator-inhibitor system which 
have a very rich behaviour. Other examples of these systems are given below. 

Example 1 

This model qualitatively incorporates activation and inhibition: 

— = cu, (4.41 

dt b + v ' ^ ' 

dv 

— = du-ev. (4.42) 

Example 2 

This model is commonly referred to as the Gierer-Meinhardt model [3]: 

du , , , 

— = a-bu^ , (4.43) 

at V 

— = u^-v. (4.44) 
at 
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Example 3 

This model is commonly referred to as the Thomas model [8]. Proposed in 1975, it is an 
empirical model based on a specific reaction involving uric acid and oxygen: 



— = a - u - pR{u,v), (4.45) 
dv 

— = a{b-v) - pR{u,v), (4.46) 

where 

represents the interactive uptake. 



Chapter 5 

Introduction to spatial variation 



We have initially considered biological, biochemical and ecological phenomena with neg- 
ligible spatial variation. This is, however, often not the case. Consider a biochemical 
reaction as an example. Suppose this reaction is occurring among solutes in a relatively 
large, unstirred solution. Then the dynamics of the system is not only governed by the dy- 
namics of the rate at which the biochemical react, but also by the fact there can be spatial 
variation in solute concentrations, which entails that diffusion of the reactants can occur. 
Thus modelling such a system requires taking into account both reaction and diffusion. 

We have a similar problem for population and ecological models when we wish to incor- 
porate the tendency of a species to spread into a region it has not previously populated. 
Key examples include modelling ecological invasions, where one species invades another's 
territory (as with grey and red squirrels in the UK [10]), or modelling the spread of dis- 
ease. In some, though by no means all, of these ecological and disease-spread models the 
appropriate transport mechanism is again diffusion, once more requiring that we model 
both reaction and diffusion in a spatially varying system. 

In addition, motile cells can move in response to external influences, such as chemical 
concentrations, light, mechanical stress and electric fields, among others. Of particular 
interest is modelling when motile cells respond to gradients in chemical concentrations, a 
process known as chemotaxis, and we will also consider this scenario. 

Thus, in the following chapters, we will study how to model such phenomena and how 
(when possible) to solve the resulting equations in detail, for various models motivated 
from biology, biochemistry and ecology. 

References. 

• J. D. Murray, Mathematical Biology Volume I, Chapter 11 [8]. 

• N. F. Britton, Essential Mathematical Biology, Chapter 5 [1]. 
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5.1 Derivation of the reaction-diffusion equations 

Let i S {1, . . . ,m}. Suppose the chemical species Cj, of concentration q, is undergoing a 
reaction such that, in the absence of diffusion, one has 

dcj 

= Ri{ci,C2,...,Cm)- (5.1) 

Recah that Ri{ci,C2, ■ ■ ■ ,Cm) is the total rate of production/destruction of Cj per unit 
volume, i.e. it is the rate of change of the concentration q. 

Let t denote time, and x denote the position vector of a point in space. We define 

• c{x,t) to be the concentration of (say) a chemical (typically measured in mol m~'^). 

• q{x,t) to be the flux of the same chemical (typically measured in mol m~^ s~^)- 

Recall that the flux of a chemical is defined to be such that, for a given infinitesimal 
surface element of area dS and unit normal n, the amount of chemical flowing through 
the surface element in an infinitesimal time interval, of duration dt, is given by 

n • q dSdt. (5.2) 




Definition. Pick's Law of Diffusion relates the fiux q to the gradient of c via 

q = -DVc, (5.3) 
where D, the diffusion coefficient, is independent of c and Vc. 



Bringing this together, we have, for any closed volume V (fixed in time and space), with 
bounding surface dV, 



dt 

Hence 



^[cidV = -[ q-ndS+ [ Ri{ci,C2, ■ ■ ■ ,0^) dV, iG {!,..., m}. (5.4) 
Jv JdV JV 



^^Qdy = J^V ■qdV + J^Riici,C2,...,Cm)dV 



(5.5) 



/ {V ■{DVci)+Ri{ci,C2,...,Cm)} dV, (5.6) 
Jv 
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and thus for any closed volume, V, with surface dV, one has 

^1^- V-pVQ)-i?i|dy = 0, {!,..., m}. (5.7) 

Hence 

(JC' 

= V ■{DVci) + Ri, xeV, (5.8) 

ot 

which constitutes a system of reaction-diffusion equations for the m chemical species in 
the finite domain D. Such equations must be supplemented with initial and boundary 
conditions for each of the m chemicals. 

Warning. Given, for example, that 

r-2TT 

/ cos6'd6l = 7^ cos 61 = 0, 6'G[0,27r], (5.9) 
Jo 

are you sure one can deduce equation (5.8)? 
Suppose 

(JC ' 

-^-V.{DVc,)-R,^0, (5.10) 

at some x = x*. Without loss of generality, we can assume the above expression is positive 
i.e. the left-hand side of equation (5.10) is positive. 

Then 3 e > such that 

_2 _ V . (DVq) - i?, > 0, (5.11) 

for aU X G Be{x*). 
In this case 



Br 

^ - V • (DVq) - R^ 



'B,{X*) 

contradicting our original assumption, equation (5.7). 



dV > 0, (5.12) 



Hence our initial supposition is false and equation (5.8) holds for x £ T>. 

Remark. With one species, with a constant diffusion coefficient, in the absence of reac- 
tions, we have the diffusion equation which in one dimension reduces to 

dc „ d^c 

For a given length scale, L, and diffusion coefficient, D, the timescale of the system is T = 
L'^/D. For a ceh, L ~ 10~^m = 10"^ cm, and for a typical protein D ~ 10 cm^s would 
not be unreasonable. Thus the timescale for diffusion to homogenise spatial gradients of 
this protein within a cell is 



L2 lO"*' cm2 

r~ — 5 ^~10s, (5.14 

D 10-7 cm2 s~i ^ ' 



therefore we can often neglect diffusion in a cell. However, as the scale doubles the time 
scale squares e.g. L x 10 ^ T x 100 and L x 100 ^ T x 10^. 
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Note. The above derivation generalises to situations more general than modelling chem- 
ical or biochemical diffusion. For example, let I{x, y, t) denote the number of infected 
people per unit area. Assume the infectives, on average, spread out via a random walk 
mechanism and interact with susceptibles, as described in Section (6.2.1). One has that 
the flux of infectives, Qj, is given by 

qj = -DiVI, (5.15) 

where Dj is a constant, with dimensions of (length)^ (time)"^. Thus, one has, via precisely 
the same ideas and arguments as above, that 

dl 

— = \/.(DjVI) + rIS-aI, (5.16) 

where S{x, y, t) is the number of susceptibles per unit area, and r and a have the same 
interpretation as in Section 6.2.1. 



Fisher's Equation. A very common example is the combination of logistic growth and 
diffusion which, in one spatial dimension, gives rise to Fisher's Equation: 

Dirj, + ru(l--], (5.17) 



dt dx^ V K 

which was first proposed to model the spread of an advantageous gene through a popula- 
tion. See Section 6.1 for more details. 



5.2 Chemotaxis 

As briefly mentioned earlier, motile cells can move in response to gradients in chemical 
concentrations, a process known as chemotaxis. This leads to slightly more complicated 
transport equations, as we shall see. 

The diffusive flux for the population density of the cells, n, is as previously: J d = —DVn. 
The flux due to chemotaxis (assuming it is an attractant rather than a repellent) is taken 
to be of the form: 

Jc = nx{a)Va = nV^{a), (5.18) 

where a is the chemical concentration and <&(a) increases monotonically with a. Clearly 
x(o) = *^'(a); the cells move in response to a gradient of the chemical in the direction in 
which the function ^{a) is increasing at the fastest rate. 

Thus the total flux is 

Jr) + Jc = -DVn + nx{a)Va. (5.19) 



Combining the transport of the motile cells, together with a term describing their repro- 
duction and/or death, plus an equation for the chemical which also diffuses and, typically, 
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is secreted and degrades leads to the following equations 

— = V-{DVn)-V-{nx{a)Va) + f{7i,a), (5.20) 
da 

— = V • (DaVa) + An - fxa. (5.21) 

In the above the above /(n, a) is often taken to be a logistic growth term while the function 
x(o) describing the chemotaxis has many forms, including 

X(a) = ^2, (5,22) 

X(a) = (5.23) 

where the latter represents a receptor law, with $(a) taking a Michaelis-Menten form [5]. 



Chapter 6 

Travelling waves 



Certain types of models can be seen to display wave-type behaviour. Here we will be 
interested in travelling waves, those that travel without change in shape and at constant 
speed. 

References. 

• J. D. Murray, Mathematical Biology Volume I, Chapter 13 [8]. 

• J. D. Murray, Mathematical Biology Volume II, Chapter 1 [9]. 

• N. F. Britton, Essential Mathematical Biology, Chapter 5 [1]. 

6.1 Fisher's equation: an investigation 

Fisher's equation, after suitable non-dimensionalisation, is 

^ = ^ + (6.1) 

where (3, z, t are all non-dimensionalised variables. 

Clearly the solution of these equations will depend on the initial and boundary conditions 
we impose. We state these conditions for the time being as 

/3(2;, t) — )• /3±oo as z — )• ±00 and f3{z,T = 0) = /3q{z), (6.2) 

where /3±oo, Po, are constants. 

6.1.1 Key points 

• We will investigate whether such a wave solution exists for the above equations which 
propagates without a change in shape and at a constant (but as yet unknown) speed 
V. Such wave solutions are defined to be travelling wave solutions. 
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The investigation of the potential existence of a travelhng wave solution will be 
substantially easier to investigate on performing the transformation to the moving 
coordinate frame y = z — vt as, by the definition of a travelling wave, the wave 
profile will be independent of time in a frame moving at speed v. 

Using the chain rule and noting that we seek a solution that is time independent 
with respect to the y variable, we have 

d£_d£dy d£dT dp_ _ d^dy_ d£dT_ 

dt dy dt dr dt dz dy dz dr dz' 

i.e. 

dt dy dr dz dy' 

Assuming /3 = I3{y) so that d/S/dr = the partial differential equation, (6.1), reduces 
to 

/3" + v(3' + /3(1 - /?) = where ' = ^. (6.5) 

dy 

One must choose appropriate boundary conditions at ±00 for the travelling wave 
equations. These are the same as the boundary conditions for the full partial differ- 
ential equation (but rewritten in terms of y), i.e. 

P{y) P±oo as y ^ ±00, (6.6) 
where /3±oo) are the same constants as specified in equation (6.2). 
One must have that /3+cxd, /3-oo only take the values zero or unity: 

[/3" + vP' + /3(1 - /?)] dy = 0, (6.7) 

J —00 
gives 

/oo 
/3(l-/3)dy = 0. (6.8) 
-00 

If we want /3 — )■ constant as y — )■ ±00 and /3, /3' finite for Vy we must have either 
/3— T-Oor/?— T-lasy— )-oo and similarly for y — )■ —00. 

• With the boundary conditions (/?(— 00), /3(oo)) = (1,0), we physically anticipate 
u > 0. 




Indeed there are no solutions of the Fisher travelling wave equations for these bound- 
ary conditions and f < 0. 



Chapter 6. Travelling waves 



50 



• Solutions to equations (6.1) and (6.2) are unique. The proof would be an exercise 
in the theory of partial differential equations. 

• The solutions of the travelling wave equations are not unique. One may have solu- 
tions for different values of the unknown v. Also, if f3(y) solves (6.5) for any fixed 
value of V then, for the same value of v, so does /?(?/ + A), where A is any constant. 
For both v and A fixed the solution of the travelling wave equations are normally 
unique. 

• Note that the solutions of the travelling wave equations, (6.5), can only possibly 
be solutions of the full partial differential equation, when considered on an infi- 
nite domain. [Realistically one requires that the length scale of variation of the 
system in question is much less than the length scale of the physical domain for 
a travelling wave to (have the potential to) be an excellent approximation to the 
reaction-diffusion wave solutions on the physical, i.e. finite, domain]. 

• One "loses" the partial differential equation initial conditions associated with (6.1) 
and (6.2). The solution of the travelling wave equations given above for P are only 
a solution of the full partial differential equation, (6.1), for all time if the travelling 
wave solution is consistent with the initial conditions specified in (6.2). 

• However, often (or rather usually!), one finds that for a particular choice of v the 
solutions of the full partial differential equation system, (6.1) and (6.2), tend, as 
t — )• CO, to a solution of the travelling wave equations (6.5), with fixed v and A, for 
a very large class of initial conditions. 

• The Russian mathematician Kolmogorov proved that solutions of the full partial 
differential equation system, (6.1) and (6.2), do indeed tend, as t — )■ oo, to a solution 
of the travelling wave equations for v = 2 for a large class of initial conditions. 

6.1.2 Existence and the phase plane 

We will investigate the existence of solutions of Fisher's equation, equation (6.5), with the 
boundary conditions (/3(— oo), /3(oo)) = (1, 0) and u > 0, by means of an extended exercise 
involving the phase plane (/?', /3). 

Consider the travelling wave equation 

d^/3 d/3 

^+.^ + /3(l-/3) = 0, (6.9) 
with V > and the boundary conditions (/?(— oo), /3(oo)) = (1,0). 

Exercise 1. Show that the stationary point at (/?,/?') = (1,0) is always a saddle point 
and the stationary point at (/3, /?') = (0, 0) is a stable node for v > 2 and a stable focus 
for V < 2. 
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Solution. Writing /3' = 7 gives 



dy I 7 / dy\ (3' [ -vj- /3(1 - /3) 



(6.10) 



The Jacobian, J, is given by 



df 


dl 


f 


dj 


dg 


dg_ 


dp 


d-f 



At (0, 0) we have 



det(J - AI) = det I ^ , | \^ + vX + l = 0, (6.12) 

\ -1 -v-X 



and hence 



A = 1 . (6.13) 



Therefore: 

• if f < 2 we have A = —v/2 it ifi and hence a stable spiral; 

• if u > 2 we have A = —v/2 it /_f and hence a stable node; 

• if f = 2 we have A = — 1 and hence a stable node. 

At (1,0) we have 

det( J - AI) = det ( ^ \ A^ + uA - 1 = 0, (6.14) 

and hence 

A = (6.15) 

Therefore (1,0) is a saddle point. 

Exercise 2. Explain why solutions of Fisher's travelling wave equations must tend to 
phase plane stationary points as y — )• ±00. Hence explain why solutions of (6.9) with 
V < 2 are unphysical. 

Solution. {f3,-f) = (/3, f3') will change as y increases, unless at a stationary point. Therefore 
they will keep moving along a phase space trajectory as y — )• 00 unless the y — )■ 00 limit 
evolves to a stationary point. 

To satisfy lim^^o /3(y) = 0, we need to be on a phase space trajectory which "stops" at 
/3 = 0. Therefore we must be on a phase space trajectory which tends to a stationary 
point with /3 = as y — )• 00. 

Hence we must tend to (0, 0) as y — )■ 00 to satisfy liniy^oo /3(y) = as y — )• 00. 
An analogous argument holds as y — )• — 00. 

If v < 2 then /9 < at some point on the trajectory which is unphysical: 
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7 




13 



Exercise 3. Show that the gradient of the unstable manifold at {f3,l3') = (1,0) is given 

by ^ ^ ^ 

(6.16) 



Sketch the qualitative form of the phase plane trajectories near to the stationary points 
for V > 2. 



Solution. We require the eigenvectors of the Jacobian at (1,0): 

' 1 \ 



1 

1 -V 



1 

Q± 



Ah 



q± 



q± = X± and 1 — vq-i- = X±q±. (6.17) 



Hence 



v± 



1 



(6.18) 




Exercise 4. Explain why any physically relevant phase plane trajectory must leave 
{/3,/3') = (1,0) on the unstable manifold pointing in the direction of decreasing /3. 



Solution. Recall that, close to the stationary point, 



(3 
7 



/3* 
7* 



a-e^-yv_+a+e^+yv+. 



(6.19) 



The solution moves away from the saddle along the unstable manifold, which corresponds 
to a_. 
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def 

Exercise 5. Consider v > 2. With 7 = /5' show that for /3 € (0, 1] the phase plane 
trajectories, with gradient 

d/3 7 ^ 

satisfy the constraint dj/dfS < —1 whenever 7 = — /?. 

Solution. 



.20) 



d7 
d^ 



-v + il-/3)< {-V + 2) - (1 + /3) < -1. 



.21) 



/9=-7 



Exercise 6. Hence show that with v >2 the unstable manifold leaving (/3,/3') = (1,0) 
and entering the region /?' < 0, /3 < 1 enters, and can never leave, the region 



n=^{{P,j) I 7<0, /3g [0,1], 7>/3}. 



(6.22) 



Solution. Along £1 = {(/3,7) | 7 = 0, /3 G (0, 1)} the trajectories point vertically into TZ 
as 

00 as we approach £1 and 7' = -/3(1 - /3) < 0. (6.23) 



d7 



d/3 



Along £2 = {(/3,7) I ,5 = 1, 7 E (-1,0)} we have 



d7 
dp 



C2 



/3(l-/3) 
7 



-V < 0. 



(6.24) 



Hence trajectories that enter TZ cannot leave. There any trajectory must end at a station- 
ary point and trajectories are forced to the point (/3,7) = (0,0). 



Exercise 7. Thus prove that that there exists a monotonic solution, with /? > 0, to 
equation (6.9) for every value of > 2 and, with v > 2 fixed, the phase space trajectory 
is unique. 

Solution. The above analysis is valid for v > 2. For v fixed a trajectory enters the region 
TZ along the unstable manifold (only one unstable manifold enters TZ). The solution is 
monotonic as 7 < throughout TZ. 

Figure 6.1 shows the results of numerical simulation of the Fisher equation (6.1) with 
initial and boundary conditions given by (6.2) at a series of time points. 



6.1.3 Relation between the travelling wave speed and initial conditions 

We have seen that, for v fixed, the phase space trajectory of Fisher's travelling wave 
equation is unique. The non-uniqueness associated with the fact that if /3(y) solves Fisher's 
travelling wave equation then so does /3{y + A) for A constant simply corresponds to a 
shift along the phase space trajectory. This, in turn, corresponds simply to translation of 
the travelling wave. 
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Figure 6.1: Solution of the Fisher equation (6.f) witli initial and boundary conditions given 
by (6.2) at times t = 10, 20, 30, 40, 50. 

Key question. Do solutions of the full system, equations (6.1) and (6.2), actually evolve 
to a travelling wave solution, and if so, what is its speed? 

Non-Examinable: initial conditions of compact support 

Kolmogorov considered the equation 

= g + (6.25) 

with the boundary conditions 

■0(z,t)— 7-1 as z — 7- — CX3 and ^/'(z,r)— t-O as z — t- oo, (6.26) 

and non-negative initial conditions satisfying the following: there is a K, with < K < oo, 
such that 

il;(z,T = 0) = for z>K and ilj{z,T = 0) = 1 for z < -K. (6.27) 

He proved that ^l^{z, r) tends to a Fisher travelling wave solution with f = 2 as t — t- oo. 

This can be applied to equations (6.1) and (6.2) providing the initial conditions are non- 
negative and the initial condition for /3 satisfies the above constraint, i.e. there is a K, 
with < i^T < oo, such that 

P{z, r = 0) = for z> K and f3{z, r = 0) = 1 for z < -K. (6.28) 

Under such constraints (3 also tends to a Fisher travelling wave solution with v = 2. 

6.2 Models of epidemics 

The study of infectious diseases has a long history and there are numerous detailed models 
of a variety of epidemics and epizootics (i.e. animal epidemics). We can only possibly 
scratch the surface. In the following, we consider a simple, framework model but even this 
is capable of highlighting general comments about epidemics and, in fact, approximately 
describes some specific epidemics. 
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6.2.1 The SIR model 



Consider a disease for which the population can be placed into three compartments: 

• the susceptible compartment, S, who can catch the disease; 

• the infective compartment, I, who have and transmit the disease; 

• the removed compartment, R, who have been isolated, or who have recovered and 
are immune to the disease, or have died due to the disease during the course of the 



Assumptions 

• The epidemic is of short duration course so that the population is constant (counting 
those who have died due to the disease during the course of the epidemic) . 

• The disease has a negligible incubation period. 

• If a person contracts the disease and recovers, they are immune (and hence remain 
in the removed compartment). 

• The numbers involved are sufficiently large to justify a continuum approximation. 

• The 'dynamics' of the disease can be described by applying the Law of Mass Action 



The model 

Then the equations describing the time evolution of numbers in the susceptible, infective 
and removed compartments are given by 



epidemic. 



to: 




r 



a 



(6.29) 



d5 

d^ 
d/ 

di 
dR 



rIS — al 



al. 



rIS, 



(6.30) 




(6.32) 



subject to 



S{t = 0) = 5o, I{t 



0) = lo, R{t = 0) = 0. 




Note that 






Key questions in an epidemic situation are, given r, a, 5*0 and Iq 
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1. Will the disease spread, i.e. will the number of infectives increase, at least in the 
short-term? 

Solution. 

— = —rIS =^ is decreasing and therefore S < Sq. (6.35) 

^ = I{rS -a) < I{rSo - a). (6.36) 

Therefore, if Sq < a/r the infectives never increase, at least initially. 

2. If the disease spreads, what will be the maximum number of infectives at any given 
time? 

Solution. 

d/ {rS-a) p def a 

- = ^ = -1 + ^ where p = -. (6.37) 

Integrating gives 

I + S — plnS = Iq + So — plnSo, (6.38) 

and so, noting that dl/dS = for S = p, the maximum number of infectives is given 
by 

Ima. = I ^° ^° - ^ . (6.39) 

\^ Io + So- plnSo- plnp- p Sq > p 

3. How many people in total catch the disease? 

Solution. From 2, / — t- as t — t- oo. Therefore the total number who catch the 
disease is 

R{oo) = No- S{oo) - /(oo) = No- 5(oo), (6.40) 
where S'(oo) < So is the root of 

^oo - yolnS'oo = iVo - yolnS'o, (6-41) 

obtained by setting S = Soo and No = lo + So in equation (6.38). 



6.2.2 An SIR model with spatial heterogeneity 

We consider an application to fox rabies. We will make the same assumptions as for the 
standard SIR model, plus: 

• healthy, i.e. susceptible, foxes are territorial and, on average, do not move from their 
territories; 

• rabid, i.e. infective, foxes undergo behavioural changes and migrate randomly, with 
an effective, constant, diffusion coefficient D; 

• rabies is fatal, so that infected foxes do not return to the susceptible compartment 
but die, and hence the removed compartment does not migrate. 



Chapter 6. Travelling waves 



57 



















t 






X. 








X 



























20 40 60 80 100 

S 

Figure 6.2: Numerical solution of the SIR model, equations (6.30)-(6.32), where the solid 
lines indicate the phase trajectories and the dashed line S' + / = 5o + /o- Parameters are as 
follows: r = 0.01 and a = 0.25. 



Taking into account rabid foxes' random motion, the SIR equations become 

— = -rIS, (6.42) 

_ = DV^I + rIS -al, (6.43) 

f = («-^^> 

The / and 5 equations decouple, and we consider these in more detail. We assume 
a one-dimensional spatial domain x E (—00,00) and apply the following scalings/non- 
dimensionalisations , 



h = ^, 5"* = -^, x^ = J^x, U=rSot, A = -^, (6.45) 
60 60 V rSo rSo 

where Sq is the population density in the absence of rabies, to obtain 

^ = V^I + I{S-X), (6.47) 
where asterisks have been dropped for convenience in the final expression. 

Travelling waves 

We seek travelling wave solutions with 

S{x,t) = S{y), I{x,t) = I{y), y = x - ct, c> 0, (6.48) 
which results in the system 

= cS'-IS, (6.49) 

= I" + cl' + I{S - X), (6.50) 
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where ' = d/dy. 

We assume A = a/(r5o) < 1 below. This is equivalent to the condition for disease spread 
in the earlier SIR model. 

Boundary conditions 

We assume a healthy population as ?/ — )■ oo: 

5^1 and / ^ 0, (6.51) 

and as y — )■ — oo we require 

1^0. (6.52) 

Bound on travelling wave speed 

We write S = 1 — P and linearise about the wavefront: 

-cP'-I = and I" + cl' + X). (6.53) 

The / equation decouples and analysis of this equation gives a stable focus at (/, /') = (0, 0) 
if the eigenvalues 

,= Z^±VIIWER, (0.54) 

are complex. This requires 

c > 2Vl - A. (6.55) 

Severity of epidemic 

S{oo) is a measure of the severity of the epidemic. We have / = cS' / S and therefore 

^(/' + cl) + cS' (^^^) = 0. (6.56) 

Therefore 

(/' + cl) + c(5 - A In S) = constant = c, (6.57) 
by evaluating the equation as — )• oo. 

In this case 

S'(-oo) - AlnS'(-oo) = 1, where S'(-oo) < 1, (6.58) 
gives the severity of the epidemic. 

Further comments on travelling wave speed 

Typically, the wave evolves to have minimum wave speed: 

C ^ Cmin- (6.59) 



Chapter 7 

Pattern formation 



Examples of the importance of spatial pattern and structure can be seen just about every- 
where in the natural world. Here we will be concerned with building and analysing models 
which can generate patterns; understanding how self-organising principles may lead to the 
generation of shape and form. 

References. 

• J. D. Murray, Mathematical Biology Volume II, Chapter 2 and Chapter 3 [9]. 

• N. F. Britton, Essential Mathematical Biology, Chapter 7 [1]. 

7.1 Minimum domains for spatial structure 

Consider the non-dimensionalised, one dimensional, budworm model, but with a diffusive 
spatial structure: 



We also suppose that exterior to our domain conditions are very hostile to budworm so 
that we have the boundary conditions 



where L > is the size of the domain. Note that /'(O) = r > 0. 

Question. Clearly u = is a solution. However, if we start with a small initial distri- 
bution of budworm, will we end up with no budworm, or an outbreak of budworm? In 
particular, how does this depend on the domain size? 

Solution. For initial conditions with < u{x,t = 0) ^ 1, sufficiently small, we can 
approximate f{u) by f'{0)u at least while u{x,t) remains small. Thus our equations are, 
approximately. 









(7.3) 
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We look for a solution of the form (invoking completeness of Fourier series) : 

oo 

■r-^ . /nvra 

u[x, t) = an{t) sm 1^— 



(7.4) 



n=l 



This gives that the time-dependent coefficients satisfy 

dan Dn^vr^ 



dt 



L2 



-«n + /'(0)an = O-nfln, 



and hence 



^{x, t) = ^ «n ^ exp 



n=l 



sm 



nvrx 



(7.5) 
(7.6) 



For the solution to decay to zero, we require that all Fourier modes decay to zero as t — )■ oo, 
and hence we require that 



a„ < Vn 



/'(0)-^^^<0 Vn, 



and thus that 



/'(O) < 



L < 



L2 
Dtt^ 

7W) 



def 



L 



crit- 



(7.7) 
(7.8) 



Hence there is a critical lengthscale, Lcrit-, beyond which an outburst of budworm is possible 
in a spatially distributed system. 

7.1.1 Domain size 

On first inspection one probably should be surprised to see that Lent increases linearly 
with the diffusion coefficient, i.e. diffusion is destabilising the zero steady state. 

We can further investigate how the nature of a steady state pattern depends on the 
diffusion coefficient. Suppose L > Lcrit and that the steady state pattern is of the form: 





u , 


i 

^max 












-L/2 


X - 


= 


L/2 X 



We therefore have 

= Du^^r. + f{u)- 
Multiplying by Ux and integrating with respect to x, we have 



= y DuxUxx dx + J Uxf{u) dx. 



(7.9) 



(7.10) 
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Thus we have 



-Dul + F{u) 



constant = F{umax) where F'{u) = f{u) 



(7.11) 



We can therefore find a relation between L, D, integrals of 

F{u)=^ r f{y)dy, 



and max(u), the size of the outbreak, as follows: 



2 \ 2 



D 



(7.12) 



\/ F{umax) — F{u) since x > and therefore Ux < 0. (7-13) 



Integrating, gives 



rL/2 ^ n 

2 dx = -(2L»)2 / 

Jo Ju 

L = {2D) 



and hence 



1 



-tax '\J F (Ufyiax) 

1 



\/F{Umax) - F{u) 



F{u) 



du. 



: dtt, 



Therefore Umax is a function of L/^2D and the root of equation (7.15). 




(7.14) 



(7.15) 



Figure 7.1: Numerical simulation of the u,n-L space, equation (7.15) with r = 0.6, q ~ 6.2 
and D = 0.1. 



7.2 Diffusion- driven instability 

Consider a two component system 



ut = DuV^u + f{u,v), (7.16) 
vt = D^V^v + g{u,v), (7.17) 
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for a; e 0, t G [0, oo) and O bounded. 
The initial conditions are 

u{x,0) = uq{x), v{x,0) = vo{x), (7.18) 
and the boundary conditions are either Diriclilet, i.e. 

u = ub, V = vb, X G dil, (7.19) 
or homogeneous Neumann, i.e. 

n • Vn = 0, n • = 0, for x e dU, (7.20) 
where n is the outward pointing normal on dU. 

Definition. Patterns are stable, time-independent, spatially heterogeneous solutions 
of equations (7.16)-(7.17). 



Definition. A diffusion- driven instability, also referred to as a Turing instability, 
occurs when a steady state, stable in the absence of diffusion, goes unstable when 
diffusion is present. 



Remark. Diffusion-driven instabilities, in particular, can drive pattern formation in 
chemical systems and there is significant, but not necessarily conclusive, evidence that it 
can drive pattern formation in a variety of biological systems. A key point is that this 
mechanism can drive the system from close to a homogeneous steady state to a state with 
spatial pattern and structure. The fact that diffusion is responsible for this is initially quite 
surprisingly. Diffusion, in isolation, disperses a pattern; yet diffusion, when in combination 
with the kinetic terms, often can drive a system towards a state with spatial structure. 



7.2.1 Linear analysis 

We wish to understand when a diffusion-driven instability occurs. Using vector and matrix 
notation we define 

uJ-]. F(u)J '^"•"l), oJ'i' » V (7.21) 



V 



^ 9{u,v) y ' y ^ 

and write the problem with homogeneous Neumann boundary conditions as follows: 

ut = DV'^u + F{u), (7.22) 

l(:)^("o"«>^(:)^(:i"::;). <-> 

with 

n • Vu = 0, xedn, (7.24) 
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I.e. 

n.Vu = = n.Vv x G dn. (7.25) 

Let u* be such that F{u*) = 0. Implicit in this definition is the assumption that u* is a 
constant vector. 

Let w = u — u* with \w\ <^ 1. Then we have 
duo 

— — = DV^w + F{u*) + Jw + higher order terms, (7.26) 
ot 



where 



/ df 


df \ 


du 
dg 


dv 


dg 


\ du 


dv J 



(7.27) 



u=u* 



is the Jacobian of F evaluated at u = u* . Note that J is a constant matrix. 

Neglecting higher order terms in \w\, we have the equation 

wt = DV^w + Jw, n ■ Vw = 0, x e 80.. (7.28) 

This is a linear equation and so we look for a solution in the form of a linear sum of 
separable solutions. To do this, we first need to consider a general separable solution 
given by 

w{x,t) = A{t)p{x), (7.29) 
where A{t) is a scalar function of time. Substituting this into equation (7.28) yields 

1 dA 

- — p = DV^p + Jp. (7.30) 

Clearly to proceed, with p dependent on x only, we require A/A to be time independent. 
It must also be independent of a; as ^ is a function of time only. Thus A/A is constant. 

We take A = XA, where A is as yet an undetermined constant. Thus 

A = Aoexp{\t), (7.31) 

for Aq ^ constant. Hence we require that our separable solution is such that 

[Xp - Jp- UVV] = 0. (7.32) 

Suppose p satisfies the equation 

V^p + k^p = 0, n • Vp = 0, a; G (7.33) 

where A; G M. This is motivated by the fact in one-dimensional on a bounded domain, we 
have p" + k'^p = 0; the solutions are trigonometric functions which means one immediately 
has a Fourier series when writing the sum of separable solutions. 
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Then we have 

[Xp -Jp + Dk^p] = 0, (7.34) 

and thus 

[\I - J + Dk'^]p = Q, (7.35) 
with IpI not identically zero. Hence 

det [\I-J + k'^D] = 0. (7.36) 

This can be rewritten as 

y -Qu X- gv + Di,k^ J 

which gives the following quadratic in A: 

+ [{Du + D,)k^ - {fu + g,)] A + h{e) = 0, (7.38) 

where 

hik") = DuD^k^ - {D„fu + D^g,)^ + {fuQv - 9ufv)- (7.39) 

Note 1. Fixing model parameters and functions (i.e. fixing Z)„, D^, /, (7), we have an 
equation which gives A as a function of k"^. 

Note 2. Thus, for any A:^ such that equation (7.33) possesses a solution, denoted Pk{x) 
below, we can find a A = A(/i;^) and hence a general separable solution of the form 

Aoe^^^'^'Pkix). (7.40) 

The most general solution formed by the sum of separable solutions is therefore 

J]Ao(fc2)e^('=')*Pfc(^), (7.41) 

if there are countable k"^ for which equation (7.33) possesses a solution. Otherwise the 
general solution formed by the sum of separable solutions is of the form 

j Ao{k^)e^^''"^'p,.{x)de, (7.42) 
where k'^ is the integration variable. 
Unstable points 

If, for any /c^ such that equation (7.33) possesses a solution, we find Re{X{k^)) > then: 

• u* is (linearly) unstable and perturbations from the stationary state will grow; 

• while the perturbations are small, the linear analysis remains valid; thus the per- 
turbations keep growing until the linear analysis is invalid and the full non-linear 
dynamics comes into play; 
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• a small perturbation from the steady state develops into a growing spatially hetero- 
geneous solution which subsequently seeds spatially heterogeneous behaviour of the 
full non-linear model; 

• a spatially heterogeneous pattern can emerge from the system from a starting point 
which is homogeneous to a very good approximation. 

Stable points 

If, for all fc^ such that equation (7.33) possesses a solution, we find Re{\{k'^)) < then: 

• u* is (linearly) stable and perturbations from the stationary state do not grow; 

• patterning will not emerge from perturbing the homogeneous steady state solution 
u*- 

• the solution will decay back to the homogeneous solution^. 

7.3 Detailed study of the conditions for a Turing instability 

For a Turing instability we require the homogeneous steady state to be stable without 
diffusion and unstable with diffusion present. Here we analyse the requirements for each 
of these conditions to be satisfied. 

7.3.1 Stability without diffusion 

We firstly require that in the absence of diffusion the system is stable. This is equivalent 
to 

i?e(A(0)) < 0, (7.43) 

for all solutions of A(0), as setting /c^ = removes the diffusion-driven term in equation 
(7.36) and the preceding equations. 

We have that A(0) satisfies 

A(0)2 - [/„ + g,] A(0) + [UQv - fvQu] = 0. (7.44) 

Insisting that i?e(A(0) < 0) gives us the conditions 

fu + gv < (7.45) 

fu9v-fv9u > 0. (7.46) 

'^Technical point: Strictly, this conclusion requires completeness of the separable solutions. This can 
be readily shown in ID on bounded domains. (Solutions of p" + k^p = on bounded domains with 
Neumann conditions are trigonometric functions and completeness is inherited from the completeness of 
Fourier series). Even if completeness of the separable solutions is not clear, numerical simulations of the 
full equations are highly indicative and do not, for the models typically encountered, contradict the linear 
analysis results. With enough effort and neglecting any biological constraints on model parameters and 
functions, one may well be able to find Du, D^, f, g where there was such a discrepancy but that is not 
the point of biological modelling. 
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The simplest way of deducing (7.45) and (7.46) is by brute force. 
The roots of the quadratic are given by 



w„^ ifu + 9v) ± V {fu + QvY - ^{fu9v - fvOu) , 

A(Uj± = . (7.47) 

7.3.2 Instability with diffusion 

Now consider the effects of diffusion. In addition to i?e(A(0)) < 0, we are required to 
show, for diffusion-driven instabihty, that there exists k"^ such that 

Re{\{k^)) > 0, (7.48) 

so that diffusion does indeed drive an instabihty. 

We have that \{k'^) satisfies 

a2 + [{Du + D,)k' - ifu + g,)] X + hiP) = 0, (7.49) 

where 

h{k^) = DuD.k" - {D,fu + D^g,)k^ + (/„5„ - gufv), (7.50) 

and 

a = {fu + gv)-{Du + D,)k^ <0. (7.51) 
Thus Re{X{k^)) > requires that 



Re (^a ± ^a^ - 4h{k^)j > ^ h{k'^) < 0. (7.52) 

Hence we must find k'^ such that 

h{k^) = DuD.k'' - {D,fu + D^g,)k^ + if^g, - g^h) < 0, (7.53) 

so that we have /c^ G [A;^,A;^] where hi^kj.) = 0. Figure 7.2 shows a plot of a caricature 
/i(A;2). 

This gives us that we have an instability whenever 

'a - _ B A + _ B 



k^G 



[kl,kl], (7.54) 



where 

A = Dyfu + Dugv and B = ADuD^{fugv - gufv) > 0, (7.55) 
and there exists a solution of the following 

+ k"^? = 0, n • Vp = 0, x£dn, (7.56) 

for A;2 in the above range. 
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Figure 7.2: A plot of a caricature h{k'^). 

Insisting that k is real and non-zero (we have considered the k = case above) we have 

A>0 and A"^ - B > 0, (7.57) 
which gives us that when Re{X{k'^)) > 0, the following conditions hold: 

^ > : D^fu + DuQv > 0, (7.58) 

A^-B>0: {DJu + A.5t.) > 2^DMfu9v - fv9u)- (7.59) 
7.3.3 Summary 

We have found that a diffusion-driven instability can occur when conditions (7.45), (7.46), 
(7.58), (7.59) hold whereupon the separable solutions, with k'^ within the range (7.54) and 
for which there is a solution to equation (7.33), will drive the instability. 

Key point 1. Note that constraints (7.45) and (7.58) immediately gives us that ^ 
D2. Thus one cannot have a diffusion-driven instability with identical diffusion coefficients. 

Key point 2. Prom constraints (7.45), (7.46), (7.58) the signs of /„, must be such 
that J takes the form 




(7.60) 



Key point 3. A Turing instability typically occurs via long-range inhibition, short-range 
activation. In more detail, suppose 



(7.61) 
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Then we have fu>0 and ^r^, < by the signs of J. In this case Dyf^ + D^g^ > D2 > 
Du- Hence the activator has a lower diffusion coefficient and spreads less quickly than the 
inhibitor. 

7.3.4 The threshold of a Turing instabiHty. 

The threshold is defined such that equation (7.39), i.e. 

D^Dykt - {D,U + Dugv)kl + UnQv - 9ufv) = 0, (7.62) 
has a single root, k^. 

Thus we additionally require 

A^ = B i.e. {Dyfu + Dugv? = ADMfu9v-gufv)>Q, (7.63) 

whereupon 

,2 ^ ^ ^ ^-"f^ + ^^9v ^ .X 

Strictly one also requires that a solution exists for 

V^p + /c^p = 0, n • Vp = 0, X e dn, (7.65) 
when k"^ = k^. However, the above value of k1 is typically an excellent approximation. 

7.4 Extended example 1 

Consider the one-dimensional case: 

ut = DuUxx + f{u,v), (7.66) 

vt = DyVxx + g{u,v), (7.67) 

for a; € [0, L], t G [0, 00) and zero flux boundary conditions at x = and x = L. 
The analogue of 

V^p + k'^p = 0, n • Vp = 0, xe dn, (7.68) 

is 

Pxx + k^P = 0, p'{0) = p'{L) = 0, (7.69) 

which gives us that 

TITT 

Pk{x) = Ak cos (kx) , k = —, n G {1,2, . . .}, (7.70) 
where A^ is fc-dependent in general but independent of t and x. 
Thus the separable solution is of the form 

cos (kx) , (7.71) 

k 

where the sum is over the allowed values of k i.e. 

fe = ^, nE{l,2,...}. (7.72) 
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Figure 7.3: Numerical simulation of the Giercr-Meinhardt model for pattern formation. 

7.4.1 The influence of domain size 

If the smallest allowed value of k"^ = vr^/L^ is such that 



vr 



= — > 

then we cannot have a Turing instability. 



(7.73) 



Thus for very small domains there is no pattern formation via a Turing mechanism. 
However, if one slowly increases the size of the domain, then L increases and the above 
constraint eventually breaks down and the homogeneous steady state destabilises leading 
to spatial heterogeneity. 

This pattern formation mechanism has been observed in chemical systems. It is regularly 
hypothesised to be present in biological systems {e.g. animal coat markings, fish markings, 
the interaction of gene products at a cellular level, the formation of ecological patchiness) 
though the evidence is not conclusive at the moment. 



7.5 Extended example 2 

Consider the two-dimensional case with spatial coordinates x = (x,y)^, x G [0, a], y € 
[0, 5], and zero flux boundary conditions. We find that the allowed values of /c^ are 



k^ 
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with 



/mvrxx /"mry 

Pm,n{^) = Am,n COS \ COS (^-^ 

excluding the case where n, m are both zero. 



n,m £ {0, 1,2,...}, 



(7.74) 
(7.75) 
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Suppose the domain is long and thin, b <ti a. We may have a Turing instabihty if 



— ^ + 
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G [k'i,kl] where = 0. 



(7.76) 



For b sufficiently small, this requires n = and therefore no spatial variation in the y 
direction. 

This means we have that the seed for pattern formation predicted by the linear analysis 
is a separable solution which is "stripes" ; this typically invokes a striped pattern once the 
non-linear dynamics sets in. 

For a large rectangular domain, 5 ~ a sufficiently large, it is clear that a Turing instability 
can be initiated with n, m > 0. This means we have that the seed for pattern formation 
predicted by the linear analysis is a separable solution which is "spots". This typically 
invokes a spotted pattern once the non-linear dynamics sets in. 



10 




10 




Figure 7.4: Changes in patterning as the domain shape changes. 



Figure 7.4 shows how domain size may affect the patterns formed. On the left-hand side 
the domain is long and thin and only a striped pattern results, whilst the on the right-hand 
side the domain is large enough to admit patterning in both directions. 

Suppose we have a domain which changes its aspect ratio from rectangular to long and 
thin. Then we have the following possibilities: 
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This leads to an interesting prediction, in the context of animal coat markings, that if it 
is indeed driven by a Turing instability, then one should not expect to see an animal with 
a striped body and a spotted tail. 




Figure 7.5: Animal coat markings which are consistent with the predictions of pattern 
formation by a Turing instability. 

Common observation is consistent with such a prediction (see Figure 7.5) but one should 
not expect universal laws in the realms of biology as one does in physics (see Figure 
7.6). More generally, this analysis has applications in modelling numerous chemical and 
biochemical reactions, in vibrating plate theory, and studies of patchiness in ecology and 
modelling gene interactions. 




Figure 7.6: Animal coat markings which are inconsistent with the predictions of pattern 
formation by a Turing instability. 



Chapter 8 

Excitable systems: nerve pulses 



Many cells communicate with one another via nerve impulses, also known as action po- 
tentials. Action potentials are brief changes in the membrane potential of a cell produced 
by the flow of ionic current across the cell membrane. Such type of communication is not 
limited to neurons by make occur in other cells, for example cardiac and muscle cells. 

See http://en.wikipedia.org/wiki/Action_potential and related links for more de- 
tails. 

References. 

• J. P. Keener and J. Sneyd, Mathematical Physiology, Chapter 4 and Chapter 8 [7]. 

8.1 Background 

Here we outline the background physics required to write down a model to describe a 
nerve impulse. Firstly, we note that: 

• numerous fundamental particles, ions and molecules have an electric charge, e.g. the 
electron, e~, and the sodium ion, Na~^; 

• it is an empirical fact that total charge is conserved; 

• electric charges exert electrical forces on one another such that like charges repel and 
unlike charges attract. The electric potential, denoted V, is the potential energy of 
a unit of charge due to such forces and is measured in volts; 

• a concentration of positive particles has a large positive potential, while a concen- 
tration of negative particles has a large, but negative potential; 

• electric current is defined to be the rate of flow of electric charge, measured in Amps. 
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8.1.1 Resistance 

Ohms Law, ISV = IR, holds in most situations, where AV is the change in potential, / is 
the current flowing and R, which may depend on material properties and geometries but 
not on I nor V is the resistance. 



> 

V + AV 



R 



I 
V 



Key point. Suppose one uses a wire of low resistance to connect a region with a con- 
centration of positive charges to a region with a concentration of negative charges. The 
charges will, very quickly, flow onto/off the wire until the potential is constant and there 
is no further flow of charge. 




8.1.2 Capacitance 

A simple example of capacitor is two conducting plates, separated by an insulator, for 
example, an air gap. 



+ 



V 




Connecting a battery to the plates, as illustrated, using wires of low resistance leads to 
charge flowing onto/off the plates. It will equilibrate (very quickly!); let Qeqm denote the 
difference in the total value of the charge stored on the two plates. The capacitance of the 
plates, C, is defined to be 

C = ^>0, (8.1) 

where C is a constant, independent of V. Thus the higher the capacitance, the better the 
plates are at storing charge, for a given potential. 
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V{t) 





Suppose now the potential is a function of time, V = V{t). Charge will flow on and off the 
plates in response to time dependent changes in V{t). If the time for a significant change 
in V{t) is far longer than the time it takes for the difference in the total value of the charge 
stored on the two plates, Q, to reach its equilibrium value, Qegm, as is essentially always 
the case, one has 

Q = Qeqm. = CV{t). (8.2) 

Hence, the current, J, i.e. the rate of flow of charge on/off the plates is given by 

J = Q = CV. (8.3) 

8.2 Deducing the Fitzhugh Nagumo equations 

An axon is a part of nerve cell: 




The nerve signal along the axon is, in essence, a propagating pulse in the potential dif- 
ference across the plasma (i.e. outer) membrane of the axon. This potential difference, 
V, arises due to the preferential permeability of the axon plasma membrane which allows 
potassium and sodium ions, K+ and iVa+, to pass through the membrane at rates which 
differ between the two ions and vary with V. In the rest state, V = Vrest — — 70mV 
(millivolts); in a nerve signal pulse in V rises to a peak of ~ 15mV. It is this pulse we are 
interested in modelling. 

The geometry of the axon can be treated as a cylindrical tube. An axon is axisymmetric, 
so we have no 9 dependence in any of our models of the axon. 

8.2.1 Space-clamped axon 

A common, simplifying, experimental scenario is to space-clamp the axon, i.e. to place a 
conducting wire along the axon's axis of symmetry. 
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• Thus, by conservation of charge, the total current flowing across the axon membrane 
must be zero. 

• Note that any changes in the interior due to, for example the axon allowing 
and Na^ to pass through its membrane, will occur on a much slower timescale and 
hence one has that the interior of the space-clamped axon has no spatial variation 
in its potential difference, no current flowing along the inside of the axon, and, most 
importantly, the total current flowing through the axon membrane is zero. 

The basic model for the space- clamped axon plasma membrane potential is given by 

= total transmembrane current per unit area, 
dV 

= C— + lNa + Ik + Io + Iapplied{t), (8.4) 

where 

• Iappiied{t) is the applied current, i.e. the current injected through the axon plasma 
membrane in the experiment, which is only function of time in most experimental 
set-ups. We will take Iappiied{t) = below. 
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• cdV/dt is the capacitance current through a unit area of the membrane. Recall: 



Rate of flow on/off capacitor = C dV/dt. 

Therefore rate of flow of charge per unit area of membrane = c dV/dt where c is the 
membrane capacitance per unit area. 

• iNa, Ik are the voltage dependent Na'^ and currents. Iq is a voltage dependent 
background current. 



These currents actually take complicated forms involving numerous other variables 
which satisfy complex equations, that can be simplified, if somewhat crudely. An 
excellent account is given in T. F. Weiss, Cellular Biophysics [12]. 

The resulting equations written in terms of the non-dimensional variables v = {V — 
yrest)/\yrest\ and T = t/T where T = 6 ms, the time scale of a typical nerve pulse, 
are 

dv 

e— = Av{5 - v){v - 1) - n, (8.5) 
dr 

dn 

— = -^n + v, (8.6) 
dr 

where A, 7, e, 6 are positive parameters such that A, 7 ~ 0(1), < e <^ 5 <^ 1. 

Key point. The spatially independent behaviour of a space-clamped axon is approxi- 
mated by the above Fitzhugh Nagumo equations, (8.5)-(8.6). 

8.3 A brief look at the Fitzhugh Nagumo equations 

We have 

= Av(5 - v)(v - 1) - n, (8.7) 

dr 

^ = -7n + v, (8.8) 
where A, 7, e, 6 are positive parameters such that ^,7 ~ 0(1), < e <^ 5 <^ 1. 

8.3.1 The {n,v) phase plane 

The nullclines of equations (8.5)-(8.6) are the lines where v = Q and n = 0. A plot of the 
nullclines separates the {v,n) phase plane into four regions, as shown in Figure 8.1. 

There are several things to note about the dynamics. 

• There is one stationary point which is a stable focus. 

• Thus, with initial conditions sufficiently close to the stationary point, the system 
evolves to the stationary point in a simple manner. 
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V 

Figure 8.1: The phase plane for the Fitzhugh-Nagumo equations with the v nullcline shown 
in red and the n nullcline in green. The trajectories for two different initial perturbations 
from the steady state are shown as dashed lines. Parameters are as follows: A — 1, j = 0.5, 
S = 0.1 and e = 0.001. 




• Consider initial conditions with n ~ 0, but v increased sufficiently. The system does 
not simply relax back to the equilibrium. However, one can understand how the 
qualitative behavioural of the system by considering the phase plane. 

• We anticipate that v = (V — Vrest)/\yrest \ behaves in the manner shown in Figure 8.2 
for a sufficiently large perturbation in v. 




0.01 



0.4 0.6 

time (t) 



0.00 



-0.01 




1.5 2.0 
time (t) 



Figure 8.2: Solutions of the Fitzhugh Nagumo equations with v dynamics indicated by the 
solid line and n dynamics by the dashed line. The right-hand figure shows the oscillations 
that arise for large t. Parameters are as follows: A = 1, j = 0.5, S = Q.l and e = 0.01. 



• This is essentially a nerve pulse (although because of the space clamping all the 
nerve axon is firing at once). 



Definition. A system which, for a sufficiently large perturbation from a stationary 
point, undergoes a large change before eventually returning to the same stationary 
point is referred to as excitable. 
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8.4 Modelling the propagation of nerve signals 

In the following, we generalise the ideas we have seen for modelling the plasma membrane 
potential of an axon to scenarios where this potential can vary along the axon. 

8.4.1 The cable model 

In the model we are about to develop we make following assumptions. 

• The cell membrane is a cylindrical membrane separating two conductors of electric 
current, namely the extracellular and intracellular mediums. These are assumed to 
be homogeneous and to obey Ohm's law. 

• The model has no 9 dependence. 

• A circuit theory description of current and voltages is adequate, i.e. quasi-static 
terms of Maxwell's equations are adequate; for example, electromagnetic radiation 
effects are totally negligible. 

• Currents flow through the membrane in the radial direction only. 

• Currents flow through the extracellular medium in the axial direction only and the 
potential in the extracellular medium is a function of z only. Similarly for the 
potential in the intracellular medium. 



These assumptions are appropriate for unmyelinated nerve axons. Deriving the model 
requires considering the following variables: 

• Ie{z,t) - external current; 

• Ii{z,t) - internal current; 

• J{z, t) - total current through the membrance per unit length; 

• Jioniz,t) - total ion current through the membrance per unit area; 

• V{z,t) = Vi{z,t) — Ve{z,t) - transmembrane potential; 

• ri - internal resistance per unit length; 

• re - external resistance per unit length; 

• C - membrane capacitance per unit area. 
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Veiz) 



Ie{z)dz Ie{z + dz)dz 

*^Ve{z + dz) 



cell exterior 



-AAAA/ 

rr,dz 



2TraJ^on(z)dz 



J{z)dz • 



V^{Z) 



r.,dz 



2-KaC^dz 



cell membrane 



h{z)dz 



V,{z + dz) 



J{z + dz)dz 



Ii{z + dz)dz 



cell interior 



Consider the axial current in the extracellular medium, which has resistance rg per unit 
length. We have 

dV 

Ve{z + dz)-Ve{z) = -rJ,{z)dz rj,{z) = -^, (8.9) 

where the minus sign appears because of the convention that positive current is a flow of 
positive charges in the direction of increasing z. Hence, if Ve{z + dz) > Ve{z) then positive 
charges flow in the direction of decreasing z giving a negative current. Similarly, 

TiUz) = (8.10) 
Using conservation of current, we have 

Ie{z + (lz,t) - Ie{z,t) = J{z,t)dz = Ii{z,t) - Ii{z + dz,t), (8.11) 

which gives 

J(r ^) = - _ 

dz dz 



J{z,t) = -^ = ^. (8.12) 



Hence 



and so 



(8 13) 



|^ = (r.+re)J. (8.14) 



Putting this all together gives 

^_ d{Ii + Ie) _ d fl dVe ^ 1 dVi\ _ / re + rA d'^Ve ^ 1 d'^V ^^^^^ 
dz dz \re dz Vi dz J \ VeVi J dz"^ ri dz"^ ' 



and so 



= ~ ^ = ~ 'JTY + {'re+ri)J{z,t) 

Ti dz'^ \ Ti J dz Ti \ dz'^ J 



.16) 
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We also have that 

J{z, t) = 2na (^J^on{V, z, t) + ) , (8.17) 

and finally, therefore, 

1 a^y av 



This gives an equation relating the cell plasma membrane potential, V ^ to the currents 
across the cell plasma membrane due to the flow of ions, Jion{V, z, t). 

Note 1. Note even though, physically, there is no diffusion, we still have a parabolic 
partial differential equation, so the techniques we have previously studied are readily 
applicable. 

Note 2. From the above equation one can model cell plasma membrane potentials given 
suitable initial and boundary conditions, and a suitable expression for Jjo„(z,t). 

We use the same expression for Jion{z,t), i.e. the expression for I^a + Ik + Iq as in the 
Fitzhugh Nagumo model of a space-clamped axon. 

Thus with V = {V — Vrest)/\yrest\ and X = Kz, where K is a constant, we have 

= e'^^+Av{6-v){v-l)-n, (8.19) 

dn 

_ = -^n + v, (8.20) 
dr 

where < ^,7 ~ 0(1), < e < 5 < 1. 

Note that K has been chosen so that the coefficient infront of the Vxx term is e^. This 
means, with respect to such variables, the front of the nerve pulse is extremely sharp. 
Hence, for such a scaling to exist, the extent of the nerve pulse must be less than eL, where 
L is the length of the axon; this constraint holds true for typical parameter estimates. The 
reason for the choice of this scaling is simply mathematical convenience in a travelling wave 
analysis. 

We are interested in nerve pulses, so we take the boundary conditions to be n, — )■ as 
X — )• zboo. 

We thus again have a system of parabolic partial differential equations to solve, and we 
are particularly interested in travelling pulse solutions. This entails that a travelling wave 
analysis would be most insightful. With the travelling wave coordinate y = x — ct and 
v{y) = v{x,t), n{y) = n{x,T), we obtain 

+ ec^ + Av{5 - v){v - I) - n = 0, (8.21) 
dy dy 

dzi 

c- -fn + v = 0. (8.22) 

dy 

We have < ^ ~ 0(1), < 7~^ (5,e < 1. One can readily investigate these ordinary 
differential equations to find that the travelling wave speed is unique, giving a unique 
prediction for the speed of a nerve pulse in terms of biophysical parameters. 



Appendix A 

The phase plane 



Throughout this appendix we wih be concerned with systems of two coupled, first-order, 
autonomous, non-Unear ordinary differential equations. 

Disclaimer. This material should have been covered elsewhere (for example in your 
course on differential equations) and hence below is intended to review, rather than intro- 
duce and lecture this topic. 

We can represent solutions to the equations 

X{x,y), (A.l) 

y(x,y), (A.2) 

as trajectories (or "integral paths") in the phase plane, that is the {x^y) plane. Suppose, 
for the initial condition x{t = UniUai) = xq, y{t = tinuiai) = Vo we plot, in the {x,y) plane, 
the solution of (A.l): 



dx 

di 
dy 

dt 




We can do exactly the same for all the values of {UniUai -.x initial ■, ViniUai] ■, to build-up a 
graphical representation of the solutions to the equations (A.l) and (A.2) for many initial 
conditions. This plot is referred to as the "phase plane portrait". 
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A.l Properties of the phase plane portrait 

The gradient of the integral path through the point (2:0,2/0) is given by 



dy dy j dx 



Y{x,y) 



i^(a;o,yo) 



Key point 1. Note that if ^(xo, yo) = and X(xo, yo) 7^ then 

= 0, 



dy 
dx 



(2^0, yo) 



which corresponds to a horizontal line segments in the phase plane. 

Key point 2. If y(xo, yo) 7^ and X(xo, yo) = then 

dy 



dx 



00 as (x,?/) (xo,2/o), 



which corresponds to a vertical line segment in the phase plane. 



(A.3) 



(A.4) 



(A.5) 



Key point 3. Assuming that either X(xo,yo) / or y(xo,yo) / 0, then two path 
integral curves do not cross at the point (xo, yo)- This is because under these circumstances 
dy/dx takes a unique value, i.e. the following is not possible: 




A. 2 EquiUbrium points 

Definition. A point in the phase plane where X(xo,yo) = y(2;o,yo) = is defined 
to be an equilibrium point, or equivalently, a stationary point. 

The reason for the above definition is because if (x,y) = (xo,yo) then both dx/dt and 
dy/dt are zero, and hence (x,y) do not change as t increases; hence x{t), y{t) remain at 
(xo,yo) for all time. 
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Key point 1. Integral curves cannot cross at points whicli are not equilibrium points. 

Key point 2. If an integral path ends it must end on a stationary point. 

Key point 3. As we shall see below, equilibrium points are only approached as f — )■ oo 
or t — >■ —CO. 

However, what about the gradient of integral paths at (xo,yo)? We informally have 

, 
-r- = (A.6) 

dx ^ ' 

which is not uniquely defined — the value ultimately depends on the details of how quickly 
X{x,y) and Y{x,y) approach zero as (x, y) — )• {xo,yo), and this generally depends on the 
direction upon which {x,y) approaches {xo,yo). 

A. 2.1 Equilibrium points: further properties 

Suppose the equations (A.l) and (A. 2) have an equilibrium point at {xo,yo). Thus 
X{xo,yo) = Y{xQ,yo) = 0. To determine the behaviour of integral paths close to the 
equilibrium point we write 

X = Xo+X, y = yQ + y, (A.7) 

where it is assumed that x, y are sufficiently small to allow the approximations that we 
will make below. 

By Taylor expansion, we have 

dX dX 

X{x,y) = X{xo + x,yo + y) = X{xo,yo) + x—{xo,yo) + y—{xo,yo) + h..o.t. 
dX dX 

= x—{xo,yo)+y-^{xo,yo) + h.o.t., (A.8) 
using the fact X{xo,yo) = 0. Similarly, we have 

_dY _dY 

Y{x,y) = Y{xo + x,yo + y) = Y{xo,yo) + x—{xo,yo) +y-^{xo,yo) + h.o.t., 
dX dX 

= x—{xo,yo)+y—{xo,yo)+h..o.t. (A.9) 
ox oy 

Note that xq and yo are constant, and hence have zero time derivative. Hence, by use 
of Taylor expansions and neglecting higher orders (i.e. taking x, y sufficiently small), we 
can neglect terms of the order O and hence we can write equations (A.l) and 

(A. 2) in the form 

^ = ( 1!'°' 'i 1!'°' 'i ] - Ju where u'il(']. (A.IO) 
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Definition. The matrix 



J = 



( 



) 



(A.ll) 



is defined to be tfie Jacobian matrix at the equihbrium point (a::o,yo)- 



A. 3 Summary 

The key points thus far are as follows. 

1. We have taken the full non-linear equation system, (A.l) and (A. 2), and expanded 
about one of its (possibly many) equilibrium points taken to be located at {xo,yo), 
using Taylor expansions of X{x,y), Y{x,y). 

2. We assume that we are sufficiently close to (xo,?/o) to enable us to only consider 
linear terms of the order of {x — xq), {y — yo)- 

3. In this way, we obtain a set of two coupled, linear, autonomous ordinary differential 
equations, i.e. equation (A. 10) above, which in principle we can solve! 

4. This procedure is sometimes referred to as "a linearisation of equations (A.l) and 
(A. 2) about the point (a;o,yo)"- 

5. In virtually all cases the behaviour of the linearised system is the same as the be- 
haviour of the full non-linear equations sufficiently close to the point (xq, yo)- In this 
respect one should note that the statement immediately above can be formulated 
more rigorously and proved for all the types of stationary points except: 

• centre type equilibrium points, i.e. case [3c] below; 

• the degenerate cases where Ai = and/or A2 = 0, which are briefly mentioned in 
item 2 on page (87). These stationary points can be considered non-examinable. 

The relevant theorem is "Hartmann's theorem", as discussed further in P. Glendin- 
ning. Stability, Instability and Chaos [4]. 

6. However, one should also note that the solution of the linearised equations may 
behave substantially differently from the solutions of the full non-linear equations, 
(A.l) and (A. 2), sufficiently far from (xo,yo)- 

A. 4 Investigating solutions of the linearised equations 

We now have a set of two coupled, linear, autonomous ordinary differential equations, 
(A. 10). It is useful to look for a solution of the form 



U = UqC 



(A.12) 
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for some constant, A. Substituting this into equation (A. 10) we obtain 

Xuoe^^ = Jtioe^* i.e. (J - XI) uq = 0. (A. 13) 

For a non-zero solution, we must have uq ^ (0, 0) and hence we require 

det(J-A/) = 0, (A.14) 
where / is the 2x2 identity matrix. 

This quadratic equation has two roots for A, denoted Ai, A2, which are possibly equal and 
possibly complex; these are, of course, the eigenvalues of J evaluated at the point (a;o,yo)- 

A.4.1 Case I 

Ai, A2 real, with Ai 7^ 0, A2 / 0, Ai / A2. Without loss of generality we take A2 > Ai 
below. 

We have two distinct, real eigenvalues. Let the corresponding eigenvectors be denoted by 
ei and 62- We thus have 

Jei = Aiei, Je2 = A2e2. (A. 15) 

We seek a solution of the form 

u = ^161 + ^262. (A. 16) 

Substituting this into equation (A. 10), we find, by comparing coefficients of ei and 62, 
that 

and hence 

Ai = Ai{t = 0)e^^\ ^2 = ^2(i = 0)e^2i_ ^^.18) 



Thus we have 



^ u = Ai{t = 0)e^i*ei + A2{t = 0)e^^'e2, (A.19) 

y / 



which gives us a representation of the solution of (A. 10) for general initial conditions. This 
information is best displayed graphically, and we do so below according to the values of 
A2, Ai. 

Note. The equilibrium point i.e. {x,y) = (0,0) can only be reached either as i — ?• 00 or 
t — )• —00. 

1. Ai < A2 < 0. The phase plot of the linearised equations in the {x, y) plane looks 
like one of the two possibilities in Figure A.l. 



Definition. An equilibrium point which results in this case is called a stable 
node, with the word "stable" referring to the fact that integral paths enter the 
node, i.e. the equilibrium point at (0,0). 
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Figure A.l: Possible phase portraits of a stable node. The equilibrium point in eaeh case is 
denoted by the large dot. 



2. A2 > Ai > 0. We still have 

u = Ai{t = 0)e-^i*ei + ^2(i = 0)e^^^e2. (A.20) 

However, the direction of the arrows is reversed as the signs of Ai, A2 are changed. 
The phase plane portraits are the same as in Figure A.l except the direction of the 
arrows is reversed. 



Definition. An equilibrium point which results in this case is called an unstable 
node, with the word "unstable" referring to the fact that integral paths leave the 
node, i.e. the equilibrium point at (0,0). 

3. A2 > > Ai. Once more, we still have 

u = Ai{t = 0)e-^i*ei + A2{t = 0)e^^^e2, (A.21) 

but again the phase plane portrait is slightly different — see Figure A. 2. 



Definition. An equilibrium point which results in this case is called a saddle 
point. 



Definition. The two integral paths originating from the saddle point are some- 
times referred to as the unstable manifolds of the saddle point. Conversely, the 
integral paths tending to the saddle point are sometimes referred to as the sta- 
ble manifolds of the saddle point. This forms part of a nomenclature system 
commonly used in more advanced dynamical systems theory; see P. Glendinning, 
Stability, Instability and Chaos [4]. 
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Figure A. 2: The phase portrait of a saddle point. The equilibrium point is denoted by the 
large dot. 

A.4.2 Case II 

A2, Ai real. One, or more, of the following also holds: 

A2 = Ai, Ai = 0, A2 = 0. (A.22) 

We typically will not encounter these degenerate cases in this course. We briefly note 
that behaviour of the full equations, (A.l), can be highly nontrivial when the linearisation 
reduces to these degenerate cases. Further details of such cases can be found in P. Glendin- 
ning, Stability, Instability and Chaos [4], which is on the reading list for this course. When 
Ai,A2 = 0, Hartmann's theorem doesn't hold. 

A.4.3 Case III 

A2, Ai complex. The complex eigenvalues of a real matrix always occur in complex conju- 
gate pairs. Thus we take, without loss of generality, 

X^ = a-ib = X*2, X2 = a + ib = Xl, (A.23) 
where a, b real, b ^ 0, and * denotes the complex conjugate. 
We also have two associated complex eigenvectors ei, 62, satisfying 

Jei = Aiei, Je2 = A2e2, (A. 24) 

which are complex conjugates of each other, i.e. ei = e^. 
Using the same idea as in Case I above, we have 

u = Ai{t = 0)e^i*ei + A2{t = Oje^^tg^, (A.25) 

though now, in general, Ai{t = 0), Ai, ei, A2(t = 0), A2, 62 are complex, and hence so is 
u. 
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Restricting u to be real gives 

u = Ai{t = 0)e^i*ei + Al{t = 0)e^^*e2 = Ai{t = 0)e^i*ei + {Ai(t = 0)e^i*ei)*, (A.26) 
and this is real, as for any complex number z, we have z + z* E M. 
After some algebra this reduces to 

u = e''* [Mcos(6t) + Ksm{bt)] = e' 

where M = (Mi, M2)^ , K = {Ki,K2)^ are real, constant vectors, which can be expressed 
in terms of Ai{t = 0), A2{t = 0) and the components of the eigenvectors ei and 62- 
Equivalently, we have 

X = e"* [cos(6t)Mi + s\Ti{bt)Ki] , y = e"* [cos(6t)M2 + sm{ht)K2\ , (A.28) 

where Mi, M2, Ki, K2 are real constants. 

1. a > 0. We have x, y are, overall, increasing exponentially but are oscillating too. 
For, example, with Ki = 0, Mi = 1 we have x = e"* cos{bt), which looks like: 




2 S h"*') 



(A.27) 



Note that the overall growth of x is exponential at rate a. Thus, in general, the 
phase plane portrait looks like one of the examples shows in Figure A. 3. 



Note. The sense of the rotation, clockwise or anti-clockwise, is easily determined 
by calculating dy/dt when y = or dx/di when x = 0. 



Definition. An equilibrium point which results in the above, is called an un- 
stable spiral or, equivalently, an unstable focus. The word "unstable" refers to 
the fact that integral paths leave the equilibrium point. 
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Figure A. 3: Possible pliase portraits of a focus. Tlie equilibrium point in each case is denoted 
by the large dot. 

2. a < 0. This is the same as 1. except now the phase plane portrait arrows point 
towards the equilibrium point as x and y are exponentially decaying as time increases 
rather than exponentially growing. 



Definition. An equilibrium point which results in this case is called a stable 
spiral or, equivalently, a stable focus. The word "stable" refers to the fact that 
integral paths enter the equilibrium point. 

3. a = 0. Thus we have A2 = —ib = — Ai, 6 7^ 0, 6 real, and 

X = [cos{bt)Mi + sm{bt)Ki] , y = [cos(6t)M2 + sui{bt)K2] , (A.29) 

where Mi, M2, Ki, K2 are constants. Note that 

K2X - Kiy = Lcos{bt), - M2X + Miy = L sm{bt) , (A. 30) 

where L = K2M1 - K1M2. Letting x* = K2X - Kiy and y* = -M2X + Miy, we 
have 

(x*)2 + (y*)2 = (A.31) 

i.e. a circle in the {x*,y*) plane, enclosing the origin, which is equivalent to, in 
general, a closed ellipse, in the (x,y) plane enclosing the origin. 

Note. As with 3. above, the sense of the rotation, clockwise or anti-clockwise, is 
easily determined by calculating dy/dt when y = or dx/dt when x = 0. 



Definition. An equilibrium point which results in this case, is called a centre. 
A centre is an example of a limit cycle. 
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Definition. A limit cycle is an integral path which is closed (and which does not 
have any equilibrium points). 



A. 5 Linear stability 

Definition. An equilibrium point is linearly stable if the real parts of both eigenvalues 
Ai, A2 are negative. 

Prom the expressions for u above, for example 

u = Ai{t = 0)e^i*ei + A2{t = 0)e^''^e2, (A.32) 

when Ai, A2 real, we see that any perturbation away from the equilibrium decays back to 
the equilibrium point. 

Definition. An equilibrium point is linearly unstable if the real parts of at least one 
of the eigenvalues Ai, A2 is positive (and the other is non-zero). 

Other situations are in general governed by the non-linear behaviour of the full equations 
and we do not need to consider them here. 

A. 5.1 Technical point 

The behaviour of the linearised equations and the behaviour of the non-linear equations 
sufficiently close to the equilibrium point are guaranteed to be the same for any of the 
equilibrium points [I 1-3], [III 1,2] or [II] with Ai = A2 7^ 0. All these equilibrium points are 
such that i?e(Ai) Re{\2) 7^ 0. This is the essence of Hartmann's theorem. This guarantee 
does not hold for centres, or the equilibrium points described in [II] with A1A2 = 0. 

The underlying reasons for this are as follows. 

• First, note from the above that integral paths which meet the equilibrium point 
can either grow/decay at exponential rate i?e(Ai), or exponential rate Re{\2), or 
consist of the sum of two such terms. Second, note that in the above we took a 
Taylor expansion. Including higher order terms in this Taylor expansion can lead to 
a small correction for the rate of exponential decay towards or growth away from 
the stationary point exhibited by the integral paths. These corrections tend to zero 
as one approaches the equilibrium point. 

• Consider the centre equilibrium point, which has Re{\i) = Re{\2) = 0, and an 
exponential growth/decay of zero. If the corrections arising from the Taylor series are 
always positive, the exponential growth/decay rate of all integral paths sufficiently 
near the stationary point is always (slightly) positive. Hence these integral paths 
grow exponentially away from the stationary point. However, b is non-zero, so x and 
y are still oscillating. Hence one has the non-linear equations behave like a stable 
focus. 
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• If Re{Xi) , Re{X2) 7^ 0, then for all integral paths reaching the stationary point, the 
above mentioned corrections, sufficiently close to the equilibrium point, are negligi- 
ble, e.g. they cannot change exponential growth into exponential decay or vice-versa. 
This allows one to show that stationary points with i?e(Ai) Re{\2) 7^ are guar- 
anteed to have the same behaviour for the linearised and the non-linear equations 
sufficiently close to the equilibrium point. 

A. 6 Summary 

We will typically only encounter stationary points [I 1-3], [III 1-3]. Of these stationary 
points, all but centres exhibit the same behaviour for the linearised and the non-linear 
equations sufficiently close to the equilibrium point as plotted above and in D. W. Jordan 
and P. Smith, Mathematical Techniques [6]. 
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